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ABSTRACT 


A  simple  radiative-conductive  model  is  developed  to 
investigate  the  thermal  characteristics  of  a  small  urban 
valley-  The  behaviour  of  the  model  is  studied  by  varying 
parameters  such  as  friction  velocity,  slope  inclination,  and 
soil  diffusivity. 

The  non-linear  flat-plane  boundary-layer  equations  are 
applied  to  the  valley,  and  integrated  numerically-  Since 
wind  is  not  predicted  by  the  model,  measurements  made  in  the 
lowest  5  meters  are  used-  Stability  effects  on  thermal 
diffusivity  are  accounted  for  with  similarity  theory.  The 
flat-plane  equations  are  modified  for  use  over  a  gently 
inclined  surface.  In  the  domain  of  the  slope  wind,  a 
neutral  wind  profile  is  superimposed  onto  the  plane-wind 
profile. 

It  is  found  that  the  inclusion  of  the  neutral  slope 
layer  explains  4  of  the  5  C°  observed  slope-rim  temperature 
difference.  This  compares  to  a  2  C°  difference  predicted 
using  the  unmodified  plane  theory.  The  effect  of  decreasing 
the  friction  velocity  is  to  decrease  evening  cooling  at  the 
1  meter  level  by  as  much  as  2  C°  while  simultaneously 
increasing  the  inversion  intensity  between  the  1  meter  level 
and  the  surface.  Decreased  soil  di ff usivities  lead  to  an 
almost  1  C°  increase.  Evening  temperatures  nearly  3  C° 
lower  result  from  more  steeply  inclined  slopes. 
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CHAPTER  1 


INTRODUCTION 


1.1  Gene ral 

Many  cities  are  situated  in  a  river  valley.  The  valley 
is  often  the  heart  of  the  city,  containing  the  focus  for 
heavy  industrial,  commercial,  and  recreational  activities. 
The  air  quality  in  such  an  urban  valley  will,  among  other 
considerations,  depend  on  its  micrometeorology.  In  cities 
where  the  valley  is  used  for  parkland  and  as  a  main  transit 
corridor,  there  is  mounting  public  concern  over  the 
continued  compatibilit y  of  these  diverse  uses.  In  this 
respect,  it  is  most  important  to  identify  the  meteorological 
factors  affecting  both  the  dispersal  and  entrapment  of 
vehicle  pollutants  in  a  valley  environment. 

Under  normal  conditions,  the  concentration  of 
pollutants  is  not  greatly  influenced  by  the  presence  of 
these  valleys,  mainly  because  of  their  small  physical 
extent.  Perhaps  this  is  why  so  few  studies  have  been 
completed  on  small  urban  valleys.  However,  the  formation  of 
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an  inversion  in  the  valley  acts  to  suppress  atmospheric 
turbulence  and  thermal  mixing,  reducing  both  the  wind  speed 
and  its  variation.  This  process  allows  the  concentration  of 
pollutants  in  the  valley.  Thus  one  major  consideration  in 
determining  the  pollution  potential  of  an  urban  valley  is 
the  evolution  of  its  temperature  field.  This  thesis 
develops  a  simple  radiative-conductive  model,  and  uses  it  to 
study  temperature  in  Edmonton’s  river  valley, 

1.2  Edmon ton  1 s  River  Valley 

The  City  of  Edmonton  is  located  at  latitude  53°33'N, 
longitude  113°30'W  on  the  banks  of  the  North  Saskatchewan 
River  in  central  Alberta.  The  Edmonton  location  is  ideal 
for  conducting  a  study  of  an  urban  valley.  The  valley 
itself  is  steep  sided  and  well  defined  with  an  average  depth 
of  50  meters  and  a  width  varying  between  1000  meters  and 
1500  meters.  It  is  cut  into  a  fairly  level  uniform  plain 
that  slopes  gently  from  the  southwest  to  the  northeast 
(Figure  1.1).  Further,  it  is  removed  from  the  influence  of 
rough  terrain  and  large  water  bodies.  The  valley  is 
generally  comprised  of  parkland,  roadways,  and  residential 
areas.  Little  industrialization  is  present.  These 
particular  geographical  characteristics  favour  the 
identification  of  the  true  valley  micrometeorology. 

To  date,  thirteen  field  experiments  have  been  carried 
out  by  Hage  (1979)  in  Edmonton's  river  valley.  Two 
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Figure  1. 1  Edmonton  Topography 
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principal  sites  have  been  used  to  collect  the  data:  a 
north-south  river  section  near  Dawson  Bridge,  indicated 
as  (A),  Figure  1.1,  and  an  east-west  river  section  near  the 
High  Level  Bridge,  indicated  as  (B) ,  Figure  1.1.  In 
addition  certain  data  from  Edmonton  City  Tower,  Namao 
Airport,  Municipal  Airport,  and  International  Airport, 
indicated  as  (C)  ,  (YED) ,  (YXD) ,  and  (YEG) ,  respectively,  in 
Figure  1.  1,  are  available  to  augment  the  experimental  data. 
In  this  thesis,  attention  will  be  focussed  on  the  three 
experiments  carried  out  in  the  vicinity  of  Dawson  Bridge 
during  the  summer  of  1978.  To  ensure  that  the 
microclimatology  of  the  valley  would  dominate  the 
observations,  it  was  attempted  to  conduct  these  experiments 
on  days  with  clear  skies  and  light  winds.  Under  these 
conditions  the  boundary  layer  is  sufficiently  decoupled  from 
the  synoptic  flow  aloft. 

1.3  The  Dawson  Bridge  Site 

The  locations  of  the  experimental  observation  stations 
at  the  Dawson  Bridge  site  are  indicated  in  Figure  1.2.  Rim 
station  8  was  situated  atop  the  east  river  bank  on  an 
exposed  short  grass  location  away  from  any  obstruction. 

Slope  station  9  was  located  in  approximately  1  meter  high 
grass  of  the  more  steeply  inclined  west-facing  slope.  Slope 
station  11  was  situated  across  the  river  in  an  almost  level 
interruption  of  the  east- facing  slope.  Upslope  of  the 
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Figure  1 . 2 


The  Dawson  Bridge  Site 
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station  was  lightly  forested,  and  further  dovnslope  was  a 
small  residential  area.  Station  11  was  located  150  meters 
downslope  from  the  rim  at  a  depth  of  approximately  30  meters 
below  plain  level  (Figure  1.3).  The  average  slope  at  the 
station  was  determined  to  be  16°.  The  orientation  of  the 
valley  at  this  point  was  estimated  to  be  13°  east  of  North. 


1.4  Instrumentation 

The  instrumentation  at  each  station  is  given  in 
Table  1.1.  The  height  of  the  instrumentation  above  ground 
is  shown  in  Table  1.2.  The  thermographs  were  always  housed 
in  a  Stevenson  screen,  and  were  either  of  the  Casella1  or 
Thies2  type.  The  temperature  profiles  were  measured  with 
chromel  constantan  thermocouples  of  0.013  millimeters 
diameter3.  Winds  were  measured  with  Eimco  CSIE04  cup 
anemometers  and  Gill5  propeller  anemometers.  Mechanical 
impulse  counters6  were  used  with  the  cup  anemometers,  and 
Rustrak  chart  recorders7  with  the  propeller  anemometers.  A 
comprehensive  description  regarding  the  electronics. 


iCasella  model  T9150/9154 
2Thies  model  2.0602.00.20 

30mega  0.013  millimeter  diameter  chromel  constantan  (type  E) 
thermocouple 

4Rimco  Commonwealth  Scientific  and  Industrial  Research 
Organization  model  ASI 

5Gill  Propeller  Anemometer  models  27103  and  27004 
6German  Post  Office  message  counter  12  Volt  6  Ampere  Type 

E16. 11 

7Rustrak  Instrument  Company  model  number  214C 
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Table  1.1  Meteorological  Instrumentation 


Station  Parameter  Instrumentation 


8 

wind  speed  and 
direction 

cup  and  vane  anemometer 

screen  temperature 

thermograph 

9 

slope  wind 

propeller  anemometer 

screen  temperature 

thermograph 

1 1 

wind  speed  profile 

cup  anemometers 

slope  wind 

propeller  anemometer 

valley  wind 

propeller  anemometer 

temperature  profile 

thermocouples 

screen  temperature 

thermograph 

Table  1.2  Instrumentation  Height 


Instrumentation 

Station 

Exper ime nt 

Heigh  t 

(meters) 

t  her  mograph 

8,  9,  11 

3, 

S,  11 

1  .2 

cup 

9 

8, 

9  ,  11 

2.64 

anemometers 

1  1 

8, 

9,  11 

1.70, 

3.  84, 

5.  8  4 

propeller 

9,  11 

8# 

9,  1  1 

o 

* 

CD 

an  emometers 

thermoco  uples 

1 1 

8 

0.34  , 

2.60, 

5.0  6, 

7.10, 

9.45 

9 

2.  44, 
9.45 

4.82, 

7.10, 

11 

0.12, 

14.73 

3.  72, 

6.  89  , 
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calibration,  and  characteristics  of  these  instruments  has 
been  detailed  previously  (Paterson , 1 978) ,  and  will  not  be 
included  here,. 


1.5  Micro meteorological  Characteristics  of  the  Valley 

Under  conditions  of  light  winds  and  clear  skies, 
observations  from  different  field  experiments  conducted  in 
Edmonton's  river  valley  showed  remarkably  consistent 
results.  This  section  describes  the  general  characteristics 
of  the  valley  micrometeorology  deduced  from  these 
observations. 

Throughout  the  afternoon  lapse  conditions  prevailed. 

Air  was  well  mixed,  both  mechanically  and  thermally. 

Edmonton  City  Tower  data  indicated  that  superadiatic  lapse 
rates  (typically  2  or  3  times  the  dry  adiabatic  lapse  rate) 
extended  to  at  least  a  height  of  90  meters  (the  top  of  the 
tower) .  Temperatures  over  both  of  the  valley  slopes  were 
generally  higher  than  those  along  the  rim.  Winds  within  the 
valley  were  consistently  weaker  than  those  observed  on  the 
rim.  This  was  likely  due  to  the  sheltering  and  aspect 
characteristics  of  the  valley.  These  lighter  winds  would  be 
responsible  for  weaker  mechanical  mixing  and  less  efficient 
transport  of  heat  out  of  the  lowest  levels  of  the  boundary 
layer.  Hence  the  higher  slope  temperatures  that  were 
recorded.  Maximum  temperatures  were  recorded  in  the  late 
afternoon  or  early  evening. 


. 
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A  one-to-one  relationship,  with  considerable  scatter, 
existed  between  the  wind  speed  at  the  rim  stations  and  that 
observed  at  the  Municipal  Airport.  Wind  speed  normally 
peaked  in  the  mid-afternoon  and  then  decreased  towards 
sunset.  Wind  direction  was  not  necessarily  that  of  the 
prevailing  flow.  This  was  due  to  the  presence  of  channel 
effects  and  compensating  flows  (Paterson, 1 978) . 

By  early  evening,  the  effects  of  decreased  insolation 
became  apparent  and  the  screen  temperatures  began  declining 
within  one  to  two  hours  before  sunset.  The  air  in  the 
valley  generally  became  isothermal  prior  to  the  air  above 
the  plain.  Similarly,  after  sunset,  the  formation  of  an 
inversion  in  the  valley  preceded  that  over  the  plain. 
Further,  the  inversion  was  observed  to  be  stronger  on  the 
north-  and  east-facing  slopes.  On  some  occasions  the 
inversion  was  observed  only  on  these  slopes.  Valley 
temperatures  were  often  comparable  to  the  rural  temperatures 
reported  at  the  International  Airport,  while  rim 
temperatures  compared  favourably  with  the  Municipal  Airport 
observations  (Hage,  1972).  The  earlier  onset  of  cooling 
over  the  slope  is  presumed  to  be  caused  by  the  difference  in 
solar  radiation  extinction  between  the  slope  and  the  plain. 
After  the  sun  has  completely  set,  the  continued  maintenance 
of  the  slope-rim  temperature  difference  is  believed  to  be 
due  to  the  presence  of  the  slope  wind  which  is  responsible 
for  both  the  drainage  of  cooler  air  down  the  slope,  and  the 
presence  of  a  shallow  neutral  layer  just  above  the  slope.  A 


. 


1 1 


neutral  layer  has  a  greater  thermal  diffusivity  than  a 
stable  layer. 

Slope  winds  were  observed  on  both  walls  of  the  valley. 
The  slope  wind  began  on  the  bank  shaded  first  (the  north- 
and  east-facing  slopes)  at  or  near  the  time  of  slope  sunset. 
Speeds  of  up  to  0.9  m/s  were  recorded  for  these  winds. 

There  appears  to  be  a  weak  15  to  20  minute  periodicity 
associated  with  oscillations  in  the  recorded  speed 
(Paterson  ,  1  978)  .  Further,  it  is  possible  that  the  wind 
speed  periodicities  are  correlated  to  those  recorded  for 
temperature.  This  gives  rise  to  the  interesting  possibility 
that  there  exists  a  brake  on  the  slope  wind  due  to 
competition  between  radiational  cooling,  advection,  and 
mixing.  However,  the  resolution  of  the  temperature  data  is 
such  that  further  investigation  is  precluded.  The  slope 
wind  is  measured  at  the  0.8  meter  level,  hence  its  vertical 
extent  is  at  least  this  value.  Paterson  (1978)  has 
established  that  the  vertical  extent  can  reach  at  least 
3  meters  at  a  location  half  way  up  the  slope,  and  at  least 
15  meters  over  the  flood  plain. 

There  is  some  evidence  for  the  existence  of  a 
downvalley  wind.  Thus  there  is  the  possibility  of  a  helical 
double  vortex  circulation  below  rim  level.  If  this  is  in 
fact  true,  the  potential  for  pollutant  concentration  within 
the  valley  under  inversion  conditions  is  enhanced. 

In  addition  to  the  valley  inversion  forming  prior  to 
the  urban  inversion,  it  was  found  that  the  valley  inversion 
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was  usually  more  intense.  Inversion  gradients  up  to 
18°C/100  m  have  been  recorded  over  the  slope.  It  has  also 
been  postulated  by  Paterson  and  Hage  (1979)  that  the  urban 
inversion  can  ‘drift1  over  the  valley  giving  rise  to  a 
double  structure  inversion.  The  exact  influence  of  the 
urban  environment  on  the  micrometeorology  of  the  valley  has 
yet  to  be  ascertained. 


13 


CHAPTER  2 


THE  THEORETICAL  MODEL 


2.  1  General 

In  this  chapter  the  governing  equations  for  the 
prediction  of  temperature  by  a  mathematical  valley  model  are 
presented.  This  set  of  equations  was  solved  numerically  by 
employing  a  finite-difference  scheme  on  an  appropiate  mesh 
of  grid  points. 

The  presence  of  strong  vertical  gradients  in  the  lower 
levels  of  both  the  wind  and  temperature  profiles  is  a 
problem  when  finite  differences  are  used.  There  are  two 
approaches  available  to  overcome  this  difficulty. 

The  first,  due  to  Estoque  (1963),  is  to  introduce  a 
constant-flux  layer  into  the  first  50  meters  of  the  boundary 
layer  and  to  interpolate  the  temperature  and  wind  profiles 
below  this  height.  Interpolation  is  accomplished  by  using 
the  'minus  one-third  power  law'  (Lumley  and  Panofsky,  1964), 
or  by  using  the  flux-profile  relationships  of  Businger,  et 
al,  (1971),  or  by  some  other  appropriate  technique. 
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The  second  approach,  adopted  by  Zdunkovski  (1971)  and 
others,  is  to  employ  a  denser  network  of  grid  points  close 
to  the  ground  thereby  more  properly  accounting  for  the 
changes  occurring  near  the  surface. 

The  nature  of  the  observations  from  the  river  valley 
dictated  that  predictions  at  points  below  50  meters  would  be 
required.  Further,  model  performance  in  the  lower  levels 
was  of  greater  interest  than  in  the  upper  levels.  Hence  the 
use  of  Estoque's  constant-flux  layer  technique  with  inter¬ 
polation  in  the  region  of  interest  was  less  desirable.  Some 
authors,  for  example  Sellers  (1965)  ,  maintain  that  it  is 
only  permissible  to  assume  a  constant-flux  layer  in  the 
lowest  meter,  and  even  then  only  under  certain  conditions. 
Thus  the  second  approach  was  adopted.  In  the  first  100 
meters,  a  logarithmic  spacing  was  used.  Logarithmic  spacing 
was  employed  so  as  to  be  compatible  with  the  logarithmic 
nature  of  the  temperature  and  wind  profiles. 

The  top  of  the  boundary  layer  was  chosen  to  be  2000 
meters.  This  is  similar  to  the  value  of  2050  meters  chosen 
by  Estoque  (1963) .  Other  investigators  have  chosen 
different  values.  For  example  Zdunkowski  (1971)  chose  3000 
meters,  and  Thorpe  and  Guymer  (1977)  chose  heights  ranging 
between  600  and  1200  meters.  According  to  Shir  (1973),  a 
grid  2000  meters  high  should  exceed  the  top  of  the  real 
boundary  layer  on  most  occasions.  More  specifically,  for  a 
station  with  Edmonton's  latitude,  the  grid  will  be  adequate 
for  friction  velocities  less  than  0.5  m/s.  In  any  case,  it 
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was  assumed  that  the  height  of  the  top  of  the  boundary  layer 
chosen  would  not  be  crucial  since  interest  is  mainly 
confined  to  the  lowest  100  meters.  The  grid  spacing  between 
the  100  meter  level  and  the  top  of  the  model  is  constant  and 
equal  to  100  meters. 

2.2  The  Prediction  Equations 

For  present  purposes  the  planetary  boundary  layer  can 
be  described  using  the  Navier-Stokes  Equation,  the 
thermodynamic  energy  or  heat  equation,  the  specific-humidity 
equation,  and  the  continuity-of-mass  eguation,  in 
conjunction  with  the  equation  of  state  and  Poisson's 
eguation.  These  equations  can  be  written,  respectively,  as 


— y-  — y  —y  — y 

-  -  Vp  -  V$  - 
P 


(2.2.1) 


+  S 


(2.  2.  2) 


(2.2.3) 


V*V 


0 


(2.2.4) 
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P  =  pRT 


(2. 2.  5) 


0  =  T 


f  PQ  lR/cP 


(2.  2-6) 


where,  after  Zdunkowski  et  al.  (1976), 

Ki  is  the  parameterized  micr otu rbulent  diffusivity  for 
either  momentum,  temperature,  or  specific  humidity, 

Fn  is  the  net  flux  density  of  short-  and  long-wave 
radiation, 

I  is  the  water  vapour  flux  density  due  to  phase 
changes,  and, 

S  is  the  sum  of  heat  sources  or  sinks. 

The  remaining  symbols  have  their  usual  meanings.  It  is 
understood  that  a  repeated  Latin  index  implies  a  summation 
with  respect  to  that  index,  and  that  the  indices  take  on  the 
values  1,  2,  and  3  corresponding  to  the  familiar  x-,  y-,  and 
z-directions  in  space. 

The  necessary  kinematic  boundary  conditions  are 

(2. 2.  7) 


—y 


V(zQ)  =  0 


V(zT)  =  V 


g 


(2. 2. 8) 


where  zG  is  the  roughness  height  and  zT  is  the  height  to  the 
top  of  the  boundary  layer. 

There  are  two  methods  of  introducing  the  surface 
temperature  iEto  the  model.  The  simpler  method  is  to  use  a 
prescribed  function  for  the  surface  temperature.  This 
approach  has  been  used,  for  example,  by  Fisher  (1961).  The 
second  method,  which  has  the  advantage  of  being  more 


. 
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realistic,  is  to  compute  the  surface  temperature  from  given 
profiles  of  the  air  temperature,  0,  and  the  soil 
temperature,  Ts ,  in  conjunction  with  the  heat  balance 
equation.  The  heat  balance  equation  is  a  statement 
regarding  the  flux  densities  through  the  air-soil  interface. 
This  method  requires  the  addition  of  a  soil  layer  to  the 
model  wherein  the  soil  temperature  equation  is  given  by 


31s  =  K  ills 

9t  s  dz2 


(2.2.9) 


The  heat  balance  approach  has  been  used,  among  others,  by 
Estoque  (1963),  and  Krishna  (1968),  and  Sasamori  (1970). 

The  heat  balance  equation  is  a  consequence  of  an 
application  of  the  law  of  heat  conservation  applied  to  a 
column  of  soil.  The  column  is  bounded  at  the  top  by  the 
air-soil  interface.  The  lower  boundary  is  taken  to  be  that 
surface  below  which  vertical  heat  exchange  is  negligible. 

For  soil  it  can  be  assumed  (Sellers,  1965)  that  the  net 
horizontal  transport  of  heat  through  the  column  is  zero. 

Thus  the  net  rate  of  change  of  heat  in  the  soil  column  will 
be  equal  to  the  sum  of  the  vertical  flux  densities  across 
the  air-soil  interface.  This  amounts  to  the  insistence  of 
continuity  in  heat  flow  between  the  soil  and  air.  The 
system  is  closed  by  requiring 

Ts(z0)  =  0(zo)  (2.2.10) 

The  boundary  conditions  imposed  at  the  top  and  bottom 


- 
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of  the  system  are 


const 


const 


(2.2.  12) 


(2,  2.  11) 


where  z  ,  the  depth  to  the  soil  bottom,  is  the  level  below 

JD 

which  vertical  heat  exchange  is  small. 

2. 3  Simplification  of  the  Prediction  Equations 

The  primary  aim  of  the  model  is  to  describe  the 
evolution  of  the  temperature  field  at  and  near  the  surface. 
This  model  would  then  be  gloved  with  a  second  model 
(Stovel,  1979)  to  yield  a  time-dependent  model  capable  of 
predicting  both  the  valley  wind  and  temperature  fields. 

Thus  the  prime  concern  of  the  present  study  was  temperature 
prediction.  Therefore,  it  was  decided  to  input  observed 
winds  into  the  model  and  remove  the  momentum  and  continuity 
equations  from  the  prediction  set. 

The  potential  temperature  equation  was  simplified  with 
the  following  assumptions: 

1.  no  advection, 

2.  no  flux  divergence  near  the  ground, 

3.  no  changes  of  water  phase, 

4.  no  heat  sources  or  sinks,  and, 

5.  no  horizontal  diffusion  of  temperature. 

Therefore  the  potential  temperature  equation  becomes 


(2.3.1) 


. 
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The  model  is  mostly  concerned  with  dry  summer  days-  Thus, 
after  Sellers  (1965),  the  latent  heat  flux  density  was 
parameterized  in  terms  of  the  sensible  heat  flux  density, 
and  the  specific  humidity  equation  removed  from  the  set  of 
prediction  equations. 

Therefore  the  set  of  equations  of  direct  concern  to  the 
model  is 


li  =  <L_  f  k  1! ' 

3t  9z  h  3Z  j 


(2.3.  2) 


with 


the 


boundar  y 


=  K  15s 

s  9z2 


conditions 


(2.3.3) 


Ts(z0)  0(zo) 

(2.3.4) 

0(zT)  =  const 

(2.3.5) 

Ts  (ZSB^  =  const 

(2.3.5) 

and  the  imposition  of  the  heat  balance  eguation  to  guarantee 
continuity  of  heat  flow. 

A  further  simplification  was  made  to  the  continuity  of 
temperature  and  heat  flow  across  the  air-soil  interface. 
After  Zdunkowski  (1976),  the  height  of  the  air-soil 
interface  was  taken  as  z=0  rather  than  z=z0 . 
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2.4  The  Heat  Balance  Equation 


sw 


C02 


■wv 


SH 


■LH 


Recalling  the  previous  s imp lif ications,  the  heat 
balance  equation  takes  the  form 

Fsw  +  FC02  +  FWV  ~  FSH  "  FLH  ”  FSHS  FLW  =  °  (2-4.1) 

where  F„rT  is  the  short-wave  radiation  to  the  ground, 

is  the  long-wave  radiation  to  the  ground  due  to 
carbon  dioxide, 

is  the  long-wave  radiation  to  the  ground  due  to 
water  vapour, 

is  the  sensible  heat  flux  density  to  the 
atmosphere, 

is  the  latent  heat  flux  density  to  the 
atmosphere, 

is  the  heat  flux  density  to  the  soil,  and, 
is  the  long  wave  radiation  from  the  ground. 

The  carbon  dioxide  flux  density  can  be  expressed  as  a 
known  fraction  of  the  terrestrial  long-wave  radiation 
(Hess,  1959)  and,  for  temperature  values  near  20°C,  this 
fraction  is  approximately  0.18.  Thus 

(2.4.2) 

where  T  is  the  surface  temperature. 

The  atmospheric  long-wave  radiation  due  to  water  vapor 
can  be  evaluated  as 


SHS 


LW 


F  -  Fttt  =  -  0.82aT4 

C02  LW 


LW 


aT4  (|£)  I*  dz 
dz 


0  (2.4.3) 
where  E  is  the  emissivity  due  to  water  vapour,  and,  w  is  the 
optical  depth  due  to  water  vapour.  For  the  model  purposes. 
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the  atmospheric  long-wave  radiation  was  assumed  constant 
over  the  forecast  period,  and  integrated  graphically  using 
the  Elsasser  chart  published  by  the  U.S.  Weather  Eureau, 
After  Kondratyev  (1969),  the  atmospheric  long-wave  radiation 
intercepted  by  the  slope  as  a  fraction  of  the  hemispheric 
long-wave  radiation  is 


cos2 (a/ 2) 


where  a  is  the  inclination  of  the  slope.  The  graphical 
integration  is  described  in  Appendix  A. 

The  short-wave  flux  density  is  a  function  of  the  solar 
constant,  the  solar  hour  angle,  the  solar  declination,  the 
latitude,  orientation,  and  inclination  of  the  slope,  and  a 
combined  albedo,  turbidity,  and  diffuse  solar  radiation 
factor.  The  precise  expression  used  is 

Fsw  =  So(1_^H  sin(6) cos (£)  +  cos (6) sin(^) cos (HAN-wt-HATO)  } 


(2.4,4) 


This  expression  is  derived  in  Appendix  B.  Values  of 
1353  W/m  2  (Paltridge  and  Platt,  1976)  and  0.39 
(Houghton,  1954)  were  used,  respectively,  for  the  solar 
constant  and  for  the  combined  albedo-turbidity-diffuse  solar 
radiation  factor  for  clear  sky  radiation. 

The  expressions  for  the  sensible  heat  flux  and  the  soil 
heat  flux  densities  are  given  by 


F 


SH 


(2.4.5) 


F 


SHS 


z=o 


(2.  4. 6) 


; 
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2.5  The  Eddy  Diff usivities 

When  the  atmosphere  is  completely  at  rest,  the  transfer 
of  physical  properties  such  as  heat  and  momentum  is 
accomplished  at  the  molecular  level.  The  transfer  process 
takes  place  according  to  such  familiar  concepts  of  molecular 
theory  as  molecular  mean  free  path  and  kinematic  viscosity. 
The  major  hypothesis  in  molecular  transport  is  that  the  mean 
flux  of  a  quantity  is  directly  proportional  to  its  gradient. 
The  proportionality  in  the  case  of  momentum  transfer  is  the 
kinematic  viscosity,  and  in  the  case  of  heat  transfer  is  the 
thermometric  conductivity.  However  the  atmosphere  is  rarely 
completely  at  rest,  even  under  extremely  stable  conditions. 
When  the  atmosphere  is  turbulent,  transport  is  accomplished 
by  turbulent  transfer  processes.  The  simplest  turbulent 
transfer  theory  is  a  quasi- analogy  to  the  molecular  transfer 
process.  The  turbulent  flow  is  said  to  possess  a  granular 
structure  in  which  small  chunks  of  fluid  or  eddies  can  break 
away  from  the  mean  flow  and  act  as  vehicles  for  heat  and 
momentum  transfer  in  a  manner  similar  to  molecules  in  the 
molecular  case.  The  hypothesis  is  that  the  eddy  flux  of  a 
quantity  is  proportional  to  its  gradient.  This  theory  is 
known  as  the  Flux-Gradient  Hypothesis,  and  the  constant  of 
proportionality  as  either  the  eddy  diffusivity  or  the  eddy 
exchange  co-efficient.  The  eddy  diffusivity  for  momentum, 
also  called  the  eddy  viscosity,  is  the  analogue  of  kinematic 
viscosity.  Similarly,  the  eddy  diffusivity  for  heat,  also 
called  the  eddy  conductivity,  is  the  analogue  of 


. 
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the rmome trie  conduct vity.  The  turbulent  analogue  of  the 
molecular  mean  free  path  is  the  characteristic  eddy  length 
scale  or  mixing  length.  This  quantity  can  be  interpreted  as 
the  distance  over  which  the  significant  eddies  will  maintain 
their  identity.  The  mixing  length,  and  other  eddy 
parameters,  are  not  solely  properties  of  the  fluid  but  also 
of  the  flow.  In  other  words,  although  it  is  character ictic 
of  the  local  turbulent  mixing,  the  length  scale  is  also  a 
function  of  its  position  in  the  flow,  the  mean  velocity  of 
the  flow,  and  the  history  of  the  flow*s  turbulent  structure. 
These  charac ter istics  demonstrate  the  prime  weakness  of  the 
Flux-Gradient  Hypothesis:  turbulence  cannot  be  properly 
treated  as  a  property  of  the  fluid.  However,  according  to 
Plate  (1971),  the  theory  is  adequate  if  the  local  turbulence 
structure  develops  in  the  same  manner  as  the  mean  flow.  The 
model  proceeds  on  this  basis. 

The  diff usivities  are  generally  complicated  and  unknown 
functions.  The  usual  approach  is  to  determine  the 
diffusivity  for  momentum  under  neutral  atmospheric 
conditions  and  then  transform  it  into  any  diffusivity  under 
any  stability  condition  using  the  corresponding  universal 
profile  function.  Monin  and  Obukhov  (1954)  have  postulated 
the  existence  of  such  universal  profile  functions,  and 
several  experimental  determinations  of  these  functions  have 
been  published  by  various  authors,  for  example  Webb  (1970), 
Paulson  (1970),  Businger,  et  al.  (1971),  and  Carl  (1973). 

Several  specifications  of  the  neutral  exchange 
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coefficient  have  been  suggested.  Ekman  (1905)  proposed 
to  be  constant  and  independent  of  height,  and  derived, 
for  the  ocean  surface  layer,  the  now  well  known  'Ekman 
spiral'  velocity  distribution.  For  present  purposes  there 
exist  three  possibilities  for  the  evaluation  of  the  neutral 
exchange  coefficient. 

The  first  expression,  due  to  Blackadar  (1962),  is 


f  kz  I2 

3Ve 

1  +  kz /X  J 

- O 

3z 

where  X  =  2- 7*1 0- * J  | /f  is  the  free  atmosphere  eddy  mixing 

length,  and,  k  is  von  Karman's  constant.  Free  atmosphere 
mixing  length  refers  to  the  characteristic  size  of  an  eddy 
at  heights  where  surface  friction  is  no  longer  an  important 
consideration,  ie.  the  eddy  is  free  of  surface  constraints. 

Physically,  the  e f f ect iveness  of  turbulent  transport  is 
determined  by  two  factors:  shear  and  eddy  size.  A  larger 
eddy  will  act  over  a  greater  distance,  and  a  larger  shear 
will  make  available  to  the  eddy  more  turbulent  energy  for 
the  transfer  process.  Blackadar's  expression  for  KN 
reflects  these  two  factors  if  the  nracketed  quantity,  which 
is  dimensionally  a  length,  is  interpreted  as  the 
characteristic  eddy  size.  The  behaviour  of  the  length  scale 
is  appropriate.  Near  the  ground,  the  eddy  size  is  of  the 
same  order  of  its  distance  from  the  ground;  and,  at  greater 
heights,  the  eddy  size  approaches  a  fixed  value. 

The  second  possible  formula  for  the  eddy  viscosity. 
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based  on  the  KEYPS6  (Sellers,  1961)  expression,  is 

K  =  kzu*  {1-eRi}  ^  (2-  5-  2) 

where  uA  is  the  friction  velocity  =  (x/p )^, 
x  = surface  stress, 

R±  is  the  Richardson  number  =  /  (f^)  #  and, 

e  is  a  constant  determined  experimentally  and  ranging 
in  value  from  9  (Pandolfo,  1966)  to  18  (Panofsky, 
Blackadar,  and  McVehil,  1960), 

The  KEYPS  equation  is  really  an  interpolation  formula  for 
the  diffusivity  in  both  the  neutral  and  unstable  regimes. 

The  first  term  is  the  contribution  due  to  forced  convection 
arising  from  mechanical  mixing.  The  second  term  represents 
the  contribution  due  to  free  convection  arising  from 
buo  yancy . 

The  last  parameterization  is  due  to  Shir  (1973)  and  is 
given  by 

KN  -  u*  'f  {  e-4z/h+  (  1  +  IfiCz/h)1-6  r1  }  (2_5_3) 

where  h  =  . 455u*/f  is  the  height  of  the  boundary  layer. 

This  particular  parameterization  arose  from  one  of  the  first 
successful  attempts  to  model  atmospheric  turbulent  flow  in 
the  planetary  boundary  layer.  The  equation  represents  the 
best  curve  fit  to  the  numerically  computed  neutral  eddy 


8from  the  initial  letters  of  the  investigators  who 
independently  proposed  the  formula:  Kazanski  and 
Monin  (1956),  Ellison  (1957),  Yamamoto  (1959), 
Panofsky  (1961) ,  and  Sellers  (1962) 


- 
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dif f usivity. 

The  present  model  does  not  predict  the  winds  at  every 
level  but  uses  instead  observed  winds  at  distinct  levels. 
Hence  a  choice  of  the  least  wind-sensitive  parameterization 
for  the  exchange  coefficient  seemed  the  most  desirable. 

Thus  Shir's  parameterization  was  adopted. 


2.6  The  Stable  and  Unstable  R  eqimes 

It  is  assumed  that  the  d if f usi vi ties  can  be  extended 
from  the  neutral  regime  by 

Ki  =  %  /  4>i(C)  (2. 6.  1) 

where  ?  =  z/L  is  the  dimensionless  height, 

L  is  the  characteristic  Obukhov  length,  and, 

^  is  the  Monin-Obukhov  profile  function  for  either 


temperature  or  momentum. 

It  is  further  assumed  that  expressions  determined  by 


Businger,  et  al.  (1971) 

from  the  Kansas  data  for 

are 

applicable 

V?<o)  - 

(1  -  150“* 

(2.6.2) 

*  (;>0)  - 
m 

1  +  4.7  ? 

(2. 6.3) 

<i>h(c<o)  = 

.74(1  -  9 

(2.  6. 4) 

<J>h(C>0)  = 

.74  +  4.7? 

(2.6.5) 

These  relations  also  provide  an  expression  for  the 
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Richardson 


number ; 

Ri(?)  = 


namely 

.74  (1 

(1  - 


C(.74  +  4.7Q 

(1  +  4 . 7  C)  z 


%  <  0 


R.  >  0 

l 


(2,  6.6) 


(2.  6.7) 


In  the  model,  the  Richardson  number  is  computed,  and  then 
used  to  determine  the  non-dimensional  height  s  from  the 
previous  equations.  The  value  of  5  obtained  is  used  to 
calculate  the  Monin-Ob ukhov  functions,  and  hence  the 
diff usivities,  for  the  appropriate  stability  regime. 


2.7  Method  of  Solution 


To  recapitulate  and  clarify  the  general  method  of  the 
model,  the  following  steps  to  solution  are  presented: 

1.  using  initial  observed  temperature  profiles  for  the 
air  and  soil,  calculate  the  surface  temperature  via 
the  heat  balance  equation, 

2.  advance  the  profiles  in  time  at  all  levels  but  the 
surf ace, 

3.  calculate  the  surface  temperature  at  the  new  time 
from  the  new  profiles  for  that  time,  and, 

4-  repeat  steps  2  and  3  for  desired  length  of  forecast. 


In  the  model,  only  the  surface  is  allowed  to  absorb  and 
re-emit  radiation.  The  remainder  of  the  boundary  layer  is 
assumed  transparent  to  radiation. 


'* 

. 


28 


According  to  Estogue  (1963)  the  assumption  of 
transparency  is  reasonable  since  the  radiative  temperature 
changes  in  the  boundary  layer  are  insignificant  in 
comparison  with  the  contributions  to  temperature  change  by 
eddy  diffusion.  However,  if  the  boundary  layer  is  polluted, 
Zdunkowski  (1976)  has  shown  that  a  less  intense  and  less 
extensive  inversion  will  result  when  radiative  flux 
divergence  calculations  are  neglected. 

It  is  assumed  that  the  present  model  atmosphere  is 
unpolluted.  Hence  the  assumption  of  transparency  is 
applicable.  The  atmospheric  long-wave  flux  density  is 
therefore  deemed  to  originate  at  the  top  of  the  model 
boundary  layer,  even  though  there  might  be  some  water  vapour 
contributing  to  this  flux  density  within  the  boundary  layer 
thereby  invalidating  the  assumption  of  transparency.  Thus 
model  heating  (or  cooling)  due  to  radiation  takes  place  only 
at  the  surface  and  the  rest  of  the  model  atmosphere  is 
heated  by  ‘conduction*,  with  the  rate  of  mixing  determined 
by  the  dif f usivities  as  a  function  of  stability. 
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CHAPTER  3 


THE  COMPQTER  MODEL 


3-1  General 

In  this  chapter  the  equations  of  Chapter  2  are 
presented  in  their  numerical  form  for  use  in  computational 
algorithms.  The  model  is  one  dimensional  in  the  sense  that 
profiles  are  generated  in  the  vertical  for  a  point  on  the 
surface.  The  model  predicts  two  such  profiles:  one  for  the 
rim  of  the  valley  (plain),  and  one  for  the  slope. 

It  is  necessary  to  run  the  model  in  extended  precision. 
The  diurnal  variation  of  temperature  during  summer  at 
Edmonton  is  10  to  15  C°/day.  To  be  conser vati ve,  it  is 
assumed  that  the  smallest  rate  of  temperature  change  is  100 
times  smaller  than  this  rate.  This  corresponds  to  a  rate  of 
0.0003  C°  per  5  minute  time  step.  Again,  to  be  on  the  safe 
side,  the  iteration  scheme  assumes  that  a  zero  rate  of 
change  has  been  reached  when  the  change  is  another  order  of 
magnitude  smaller,  ie.  a  rate  of  change  of  0.00001  C°  per  5 
minute  time  step.  Potential  temperature  is  generally  near 
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300  °K.  Therefore  a  minimum  of  eight  digits  are  necessary. 
If  the  temperature  is  expressed  as  an  initial  value  plus  a 
perturbation,  a  minimum  of  seven  digits  would  still  be 
necessary.  Therefore  extended  precision  was  adopted. 


3-2  The  Grid 

As  already  indicated,  a  logarithmic  grid  was  used  below 
100  meters  and  an  egui-spaced  grid  above.  The  logarithmic 
grid  is  partitioned  into  two  distinct  sections,  one  for  the 
atmosphere  and  one  for  the  soil.  The  air-soil  interface  is 
common  to  both  sections.  The  labelling  of  the  grid  points 
is  as  in  Figure  3.1.  To  allow  flexibility,  the  grid  is 
generated  by  the  programme.  In  the  soil  layer,  it  is 
necessary  to  specify  the  depth  of  the  first  point,  the  depth 
of  the  last  point,  and  the  total  number  of  grid  points 
desired.  The  grid  points  are  then  logarithmically  spaced. 

In  the  atmospheric  layer,  it  is  necessary  to  specify 
the  height  of  a  reference  point,  the  height  of  the  top  of 
the  sub-layer,  the  height  of  the  top  of  the  boundary  layer, 
the  number  of  points  required  between  the  reference  point 
and  the  ground,  the  number  of  points  required  between  the 
reference  point  and  the  top  of  the  sub-layer,  and  the 
spacing  desired  between  points  in  the  upper  sub-layer. 

Following  Estogue  (1963) ,  the  depth  of  the  soil  layer 
was  taken  as  -0.4  meters.  The  diurnal  soil  temperature  wave 
amplitude  is  generally  small  at  this  depth.  This  fact  is 


. 
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GRID  POINT 


REFERENCE  LEVEL 


NZPTS 


NZPTSA 


4 

3 

2 

1 

2 

3 

4 


NZSPTS 


Top  of  the  Boundary  Layer 

(2000  m) 


Top  of  the  Sub-Layer 
(100  m) 


Ground  Surface 
(0  m) 


Bottom  of  the  Soil  Layer 
(-0.4  m) 


Figure  3.1  The  Model  Grid 
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borne  out  by  Lettau  and  Davidson  (1  957). 


3.3  The  Fin ite-D if fere nee  Equations 


According  to  Smith  (1965) ,  equations  of  the  form 


9¥(z,t)  =  J2^(z,t) 

3 1  9 


are  unconditionally  stable  and  convergent  in  their 
finite-difference  form 


yk+1 
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26z 
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2¥k+1  +  ^k+1  } 
3  3-1 


2¥k  +  yk  ,  } 

J  J-l 


if  a  +  b  =  1  ,  and  1  <  b  <  2.  Subscripts  and 
superscripts  are  used  to  identify  space  and  time  inc 
respectively. 

Both  the  air  and  soil  temperature  equations  are 
parabolic  partial  differential  equations  of  the  form 


9V(z,t)  _  3  f  ,  ,3V(z,t)  ^ 

3t  3z  ^  K(  }  3 z  J 

For  an  irregularly-spaced  grid,  this  equation,  in  fi 
difference  form,  becomes  (see  Figure  3.2) 
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Figure 


3.2  Finite-Differences  on  an  Irregularly  Spaced  « 
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It  is  assumed  this  finite-difference  equation  is  also 
unconditionally  stable  and  convergent  for 
a  +  b  =  1  ,  and  1  <  b  <2. 

Zdunkowski  (1976)  has  found  this  to  be  valid  for  the  values 
a  =  1.5  and  b  =  0.5.  These  numbers  were  adopted  in  this 
model. 


3.4  Calculation  of  the  Eddy  Dif f usivities 


The  neutral  exchange  coefficient  %  was  calculated 
according  to  the  finite-difference  form  of  Shir* s 
parameterization,  equation  (2.5.3),  with  the  inclusion,  for 
completeness,  of  the  molecular  viscosity.  This  expression 
is 


K- 


■Nj+h> 


 u*k 


j3+h 
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— ( 4 /h) Z 


e  '  '  +  (1  +  lecz^/h)1-6)  1  } 


+  K 


mol 


(3.4.1) 


where  h  =  0.455uv/f,  and,  K  ^  is  the  molecular  diffusivity 
of  air.  The  friction  velocity,  u*  ,  obtained  from  the 
observational  data,  was  periodically  input  into  the  model  at 
a  predetermined  time  interval  1  £*.  This  friction 
velocity  was  calculated  from  the  observations  using 
(Businger,  1973) 
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(3.4.2) 


The  model  was  started  at  local  noon.  This  was  usually 


several  hours  before  the  observed  valley  data  were 
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available.  The  earlier  model  start-up  was  desirable  for  two 
reasons.  Firstly,  if  necessary,  the  model  has  time  to 
settle  down  prior  to  the  comparison  with  the  observations. 
Secondly,  profile  initialization  can  better  be  accomplished 
at  this  time. 

The  noon  model  start-up  time  complicated  the 
calculation  of  the  friction  velocity  for  the  rim,  since 
values  for  the  mean  wind  at  2.64  meter  height  were  not  yet 
available.  To  this  end,  it  was  assumed  that  the  wind  at  the 
Municipal  Airport  could  be  taken  as  the  rim  wind  at 
2.64  meters.  The  airport  wind  is  measured  at  10  meters 
using  a  short-period  average  value  from  a  slow- response 
anemometer.  The  error  due  to  measurement  at  different 
heights  is  expected  to  be  masked  by  the  error  due  to  both 
the  removed  location  and  less  sensitive  anemometer. 

Further,  because  experience  has  shown  that  calm  conditions 
rarely  prevail  in  central  urban  districts,  and  because  of 
the  relatively  high  threshold  value  of  the  airport 
anemometer,  a  calm  reading  at  the  airport  was  assigned  the 
value  of  1  m/s.  Municipal  Airport  data  were  also  used  when 
it  was  desired  to  run  the  model  past  the  end  of  the 
observation  period.  Whenever  the  actual  rim  wind 
observations  were  available  they  were  used.  Thus  the 
friction  velocity  for  the  rim  was  calculated  by 

u*rim  0.4  u  /  log[  - ^ - )  (3.4.3) 


where  Von  Harman* s  constant  is  taken  as  0.4,  the  surface 


■ 
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roughness  as  0.01  meters,  and,  the  mean  Hind  as  either  the 
measured  valley  rim  wind  or  the  airport  wind. 

During  the  observation  periods,  comparison  between  the 
recorded  2.64  meter  wind  at  the  rim  and  the  wind 
interpolated  to  the  same  height  over  the  slope  showed  the 
slope  wind  to  be  weaker  than  the  rim  wind  on  the  average  by 
an  amount 


n 


ave {  u  n 

slope 


/  u  . 
rim 


} 


of  the  rim  wind.  It  was  assumed  that  this  fraction,  n  ,  was 
valid  outside  the  observational  time  frame.  Thus  the 
friction  velocity  for  the  slope  was  calculated  by 
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*slope 


=  0 . 4  riu 


rim 


/  £og( 


2.64  +  0.25 
0.25 


(3.  4.  4) 


where  the  slope  surface  roughness  is  taken  as  0..25  meters. 
Hence 

^2e  a  2>28n 

U*rim  (3.4.5) 

These  friction  velocities,  calculated  at  the  2.64  meter 
level,  were  used  in  calculating  the  dif f usivities,  KN^+^  , 

at  all  model  levels.  In  practice,  the  friction  velocities 
were  calculated  at  one-half  hourly  intervals  and  then 
smoothed  with  a  three  point  running  mean. 

The  kinematic  viscosity  used  for  air  (Sutton,  1953)  is 

K  ,  =  1.5*10-5  m  2/sec  (3.4.6) 

mol 

Once  the  neutral  eddy  viscosity  has  been  calculated, 
the  desired  diff usivit ies  are  determined  by  equation 
(2.6.1).  This  equation  requires  the  evaluation  of  the 
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Monin-Obukhcv  profile  functions. 


3. 5  Calculation  of  the  Monin-Cbukhov  Profile  Functions 

The  Businger  relations,  equations  (2,6,4)  and  (2.6.5), 
are  used  to  evaluate  the  Moni D-Obukho v  profile  functions  for 
heat.  The  equations  require  the  dime nsionless  height 
which  may  be  calculated,  via  equations  (2.6.6)  and  (2.6.7), 
from  the  Richardson  number.  Thus  the  Richardson  number  must 
be  computed. 

The  Richardson  number  has  previously  been  defined  as 


R-i 


(3.5.1) 


Using  this  fcrm  of  the  Richardson  number  poses  a  problem. 
Strictly  speaking  the  wind  shear  and  hence  the  wind  must  be 
known  at  every  level.  But  wind  is  not  predicted  by  the 
model.  Further,  only  observed  winds  at  one  level  are  input 
into  the  model.  To  overcome  this  difficulty,  the  bulk 
Richardson  number,  in  conjunction  with  stability  relations 
similar  to  those  of  Golder  (1972),  was  employed.  The  bulk 
Richardson  number  (Lettau  and  Davidson,  1957)  is  defined  by 


(3.5.2) 


where  z  is  the  geometric  mean  height  between  the  top  and 
bottom  of  the  layer  considered.  Only  one  wind  measurement 
is  required  to  determine  the  bulk  Richardson  number.  The 
bulk  Richardson  number  is  related  to  the  Richardson 
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number  by 
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(3.5.3) 


Also  the  Businger  profile  functions  for  momentum  may  be 
integrated  (Paulson  1970)  to  yield 
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Equations  (2.6.6),  (2.6.7),  (3.5.3),  and  (3.5.4)  may  be 

combined,  remembering  that 
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l  =  BUog(z/zQ)  -  ij;m}2 
C  (. •  74  +  4.7cJ 


C>0 


(3.  5. 6) 

These  stability  relations  differ  from  those  of  Golder  (1972) 
in  that  the  ratio  of  the  exchange  co-efficients  is  not 
assumed  to  be  unity.  There  is  only  one  unknown  in  the  above 
equations,  and  that  is  the  dimensionless  height  C  .  If  the 
bulk  Richardson  number  is  negative,  c  is  solved  for 
iteratively  using  Newton*s  Method  (Conte  and  de  Boor,  1972). 
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By  taking  the  bulk  Richardson  number  equal  to  the 
dimensionless  height  as  a  first  guess,  solution  to  within 
5  per  cent  of  r,  is  usually  encountered  in  three  iterations. 
If  the  bulk  Richardson  number  is  positive,  the  solution  for 
£  is  analytical.  With  a  value  for  C  now  determined,  the 
profile  functions,  and  hence  the  diff usivities  can  be 
calculated. 

To  guard  against  the  possibility  of  the  di ff usivities 
inappropriately  assuming  molecular  values,  a  minimum 
friction  velocity  and  hence  maximum  <t>h  is  determined  from 
the  observed  data,  and  then,  if  a  <J>h  is  computed  that 
exceeds  this  maximum,  it  is  reset  to  the  maximum  value.  In 
the  actual  computations,  the  formula  for  the  eddy  viscosity 
is  recast  (Businger,  1973)  as 

Km  =  u*  /  (  9u/3z  )  (3.5.7) 


The  exact  procedure  to  determine  the  maximum  <f>h  is  as 
follows 

1.  determine  u*min  via  equation  (2.4.2)  using  the 

minimum  u  observed, 

2.  determine  via  equation  (3.5.7), 

3.  assume  %  =  neutral  is  9iven  bY  Km  at  sunset,  and 

determine  K^  via  steps  1  and  2  using  observed 
sunset  values,  and, 

4.  using  equations  (2.6.1)  and  (2.6.3),  calculate 
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5.  hence  calculate  <f>h 


max 


■ 
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3.6  Calculation  of  the  Surface  Temperature 

The  method  used  to  calculate  the  surface  temperature 
follows  the  same  procedure  as  Jacobs  and  Brown  (1973).  The 
heat  balance  equation  is  recast  as 

0  =  T  -  {  (  Fsw  +  Fwv  -  fsh(t)  -  flh(t)  -  fshs(t)  5/-82c'  r*4 

(3.6.1) 

This  is  an  implicit  equation  in  the  surface  temperature,  T. 
Difficulties  arise  in  solving  this  equation  iteratively. 

The  quantity  in  brackets  can  take  on  values  less  than  zero. 
To  overcome  this  problem,  let  Tc  be  the  temperature  at  which 
the  quantity  in  brackets  becomes  zero,  and  define  a  function 

T  -  «FSW  +  Fwv  -  Fsh(T)  -  Flh(T)  -  FSHS(T))/.82ar\  TSTC 
T  ,  T>TC 

(3.  6.2) 

The  real  root  of  g  (T)  =  0  determines  the  surface 
temperature.  For  T  <  Tc,  the  root  is  found  iteratively 
either  by  Newton's  Method  (Conte  and  de  Boor,  1972)  or,  if 
the  root  lies  in  a  region  where  g(T)  is  not  continuously 
differentiable,  by  the  Secant  Method  (Conte  and  de 
Boor,  1972).  For  T  >  Tc,  the  Reguli  Falsi  Method  (Conte  and 
de  Boor,  1972)  is  used-  If  T  happens  to  be  greater  than 
Tc  on  the  first  iteration,  an  initial  value  of  T  =  250  °K  is 
used  in  the  Reguli  Falsi  method.  This  value  is  usually 
smaller  than  any  Edmonton  summer  surface  temperature.  Thus 
bracketing  of  the  root  is  accomplished. 


. 


Ui- 
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The  sensible  and  soil  heat  flux  densities  are  of  the 

form 


F 


z=o 


(3.6.3) 


Due  to  the  logarithmic  nature  of  temperature  near  the 
surface,  a  second  order  Taylor  Series  expansion  is  used  in 
the  finite-difference  approximation  of  the  derivative, 
namely 


3T 

9z 


z=o 


z2T(zt)  ~  Z1T(Z2)  ~  H  ~ 

Z1  z2  (z2  '  zl^ 


(3. 6.4) 


The  latent  heat  flux  density  is  parameterized  in  terms 
of  the  sensible  heat  flux  density  (after  Sellers,  1965)  in 
the  following  manner 


F 

SH* 

F 

SH* 


>  0 

<  0 


(3.6. 5 ) 


3. 7  Prediction  Procedure  for  the  Temperature  Profiles 

The  atmospheric  temperature  profile  is  approximated  by 
the  temperature  values  at  N=NZPTS  discrete  points.  It  is 
assumed  that  the  temperature  at  the  top  of  the  profile  (the 
top  of  the  boundary  layer)  is  constant.  Similarly,  the  soil 
temperature  profile  is  approximated  by  values  at  N=NZSPTS 
points,  and  the  soil  bottom  temperature  assumed  constant. 
Further,  the  temperature  at  the  interface  is  determined  by 
the  heat  balance  equation.  Thus  the  temperature  is  only 
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forecast  at  N-2  points,  namely  2  through  N-1.  Equation 
(3.  3.  4)  can  be  rewritten  as 

A.  +  B.  Vk+1  +  C.  =  D.  (3.7.1) 

3  3  + 1  J  J  3  3-1  3 

This  set  of  equations  forms  a  tridiagonal  matrix  for  the 
interior  points,  and  can  be  solved  (see  Appendix  C)  using 
Gaussian  elimination,  ie. 


yk+1 

j 


F. 

3 


G. 

3 


D. 

,1 

<1 

i 

Fi+1 

B  . 

+  A. 

G .  , 

3 

J 

3  +  1 

D 

-  A. 

F.  , 

,1 

3 

,1+1 

B 

+  A 

G.  . 

3 

J 

3+1 

-  c. 

1. 

B.  +  A.  G 
3  3  3+1 


-  C 


B  .  +  A. 
3  3 


u/k+1 

j-1 


(3.7.  2) 


(3.  7.3) 


(3. 7.4) 


Thus  a  procedure  to  find  the  profile  at  time  t  +  (k+1)At  is 

1.  assume  surface  temperature  changes  little  over  one 

k+1  lr 

time  step  and  take  T ^  =  T^, 

2.  solve  for  Tk+1  ,  j=2,...,N-1  by  equation  (3.7.2), 

j 

3.  using  the  new  profile  obtained  in  step  2,  calculate 

mk+l 

a  new  ,  and, 

4.  advance  one  time  level,  and  repeat  steps  1  to  3. 

This  procedure  has  an  inherent  weakness.  Due  to  step  1,  the 
time  step  must  be  small.  However,  this  problem  can  be 
overcome  by  employing  an  iterative  procedure.  Before 
advancing  to  step  4,  simply  repeat  steps  1  to  3  until  the 


•  £ 
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surface  temperature  calculated  at  the  most  recent  iteration 
differs  from  that  obtained  from  the  previous  iteration  by  a 
predetermined  small  amount-  As  mentioned  in  section  3.1, 
0.00001  C°  is  an  appropriate  difference,  and  this  value  was 
used.  With  this  approach,  time  steps  as  large  as  20  minutes 
have  given  consistent  profiles.  However,  the  settling-down 
time  needed  for  the  model  increases  with  the  size  of  time 
step.  A  time  step  of  5  minutes  generally  ensures  stable 
results  after  half  an  hour.  This  can  be  reduced,  if 
desired,  by  re- initializing  the  model  after  2  or  3  time 
steps.  Generally,  the  oscillations  initially  present  are 
not  severe  and  damp  out  quickly. 


3.8  Inclusion  of  a  Neutral  Slope  Layer 

By  late  evening,  the  model  consistenly  predicts  a 
one  meter  slope  temperature  that  lies  close  to  the  rim 
temperature.  This  is  contrary  to  a  3  to  6  C°  spread 
observed  when  there  is  a  slope  wind  present.  The  behaviour 
of  the  model  is  not  surprising  since  no  provision  for  the 
increased  mixing  due  to  the  slope  wind  has  been 
incorporated . 

In  Chapter  1,  it  was  indicated  that  the  slope  wind  had 
an  effective  depth  of  at  least  3  meters  and  possibly  as 
large  as  15  meters.  For  model  purposes,  the  slope  wind  was 
taken  to  be  operative  up  to  the  tenth  grid  point,  or 
approximately  4.6  meters.  It  was  then  postulated  that  the 
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mechanical  mixing,  arising  due  to  the  slope  wind,  creates  a 
neutral  layer  within  the  lowest  levels  over  the  slope  when 
it  should  otherwise  be  a  stable  layer.  To  account  for  this 
phenomenon,  the  model  was  adjusted  so  that,  at  the  moment 
the  slope  wind  started,  it  created  a  neutral  layer  below 
about  4,6  meters.  The  di ff usivities  within  this  neutral 
layer  were  then  recalculated  using  input  observed  slope 
winds.  This  procedure  was  successful  in  explaining  a  much 
larger  fraction  of  the  observed  temperature  difference. 
Presumably  advection  would  account  for  the  remainder,  but 
the  incorporation  of  an  advective  term  was  beyond  the  scope 
of  this  one  dimensional  model. 
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CHAPTER  4 


MODEL  INITIALIZATION 


4 . 1  General 

The  river  valley  experiments  carried  out  in  the  summer 
of  1978  did  not  include  measurements  of  soil  temperature  or 
measurements  of  air  temperature  above  tower  height. 

Further,  only  a  few  air  temperature  measurements  below  the 
1  meter  level  were  made.  The  model  requires  reasonably 
complete  soil  and  air  temperature  profiles  at  start  up  time. 

In  order  to  supply  this  start  up  data,  an  intelligent 
extrapolation  of  the  available  data  must  be  made.  To  this 
end,  data  from  the  Great  Plains  Turbulence  Field  Programme 
(Lettau  and  Davidson,  1957)  was  surveyed.  It  was  noticed 
that  if,  at  local  noon,  the  temperature  profiles  for  the 
first  and  third  experiment  days  were  displayed  as  the 
difference  in  temperature  from  the  temperature  at  a 
specified  reference  level,  they  represented  a  common 
log-linear  curve.  The  assumption  was  made  that  this  curve 
would  be  applicable  to  Edmonton.  Initial  guess  temperature 
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profiles  for  below  the  1.2  meter  screen  height  were  then 
extrapolated  using  screen  height  temperatures-  The  initial 
profile  between  1.2  meters  and  850  millibars  was  determined 
by  assuming  a  constant  lapse  rate  between  the  temperatures 
at  these  two  levels.  The  initial  profile  above  the  850 
millibar  level  was  determined  solely  from  the  radiosonde 
data.  It  was  assumed  that  the  850  mb  temperature  at 
radiosonde  time  would  not  differ  greatly  from  its 
temperature  at  local  noon. 

The  values  for  the  soil  thermometric  and  thermal 
conductivities  used  in  the  model  were  also  those  of  the 
Plains  data,  namely  1.5*10- 7  m2/sec  and  0.25  J/m/sec/°K. 

4 . 2  The  Atmospheric  Temperature  Pro file 

The  temperature  profile  from  the  first  experiment  of 
the  Plains  Field  Programme  is  plotted,  with  the  0.1  meter 
temperature  subtracted  from  every  value,  in  Figure  4.1. 
Superimposed  is  a  similar  plot  obtained  from  the  third 
Plains  experiment-  A  straight  line  is  fitted  visually  to 
the  data  (given  the  magnitude  of  the  assumption  made,  a  more 
accurate  fit  was  deemed  unnecessary) .  The  axis  of  the  graph 
was  then  shifted  so  that  a  temperature  difference  of  zero 
co-aligned  with  the  screen  height.  Hence  Table  4.  1  was 
constructed.  The  last  entry  in  the  table  for  0.0  meters  was 
obtained  by  linearly  extrapolating  between  0.01  and  0.05 


mete  rs. 


■ 


£°glO  (  Height  / 
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Temperature  Difference  (C°)  from 
the  Reference  Level  Temperature 


Figure  4. 1 


The  Interpolation  Curve  for  the  Initial 
Atmospheric  Temperature  Profile 
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Table  4.  1  The  Interpolation  Table  for  the  Initial 
Atmospheric  Temperature  Profile 


Height 

(m) 

Difference  from  the 

Screen  Temperature 
<C°) 

1.20 

0.0 

1.  00 

0.  2 

0.  46 

0.8 

0.  22 

1.5 

0.  10 

2.  2 

0.  05 

2.  9 

0.  00 

4.6 

■ 
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Thus  the  model  initialization  procedure  is  to  take  the 
screen  temperature  at  local  noon  and  using  Table  4. 1  arrive 
at  a  guess  for  the  temperature  profile  in  the  lowest 
1-2  meters.  The  discretization  of  the  input  profile  need 
not  align  with  the  grid.  The  model  has  the  capability  of 
interpolating  grid  point  values  from  the  input  profile.  The 
input  profile  is  converted  to  potential  temperature  by 


where 


e . 

3 


Pi  'i  •  286 
T .  — 

:  l  P3  J 


P  . 
3 


-zj  +  1)/(RT) 


(4.  2.  1) 


(4.2.2) 


T 


T  +  T 
sfc  850 


(4.  2.3) 


T  T 

_  sfc  .  +  sfc 

T  r  =  min  max 

sfc  - ^ - 


(4.2.4) 


4. 3  The  Soil  Temperature  Profile 

In  a  fashion  similar  to  that  of  Section  4.2,  the  Plains 
data  (Lettau  and  Davidson,  1957)  were  used  to  construct 
Figure  4.2  and  Table  4.2  for  soil.  In  Figure  4.2,  the  curve 
fitted  to  the  data  below  10  cm  was  justified  on  the  basis 
that  the  data  below  this  level  was  measured  by  a  different 
observer  with  different  instruments.  The  trend  of  the  curve 


is  still  maintained. 


' 


fcog10(  Depth  / 
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Temperature  Difference  (C°)  from 
the  Reference  Level  Temperature 


Figure  4.2  The  Interpolation  Curve  for  the  Initial  Soil 
Temperature  Profile 
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Table  4.2  The  Interpolation  Table  for  the  Initial  Soil 
Temperature  Profile 


Soil  Depth 

(cm) 

Difference  from  the 

Surface  Temperature 
(CO) 

0.  00 

0.  0 

0.  50 

-  2.1 

1.  04 

-  4.2 

2.15 

-  6.4 

4.  47 

-  8.7 

9.28 

-1  0.9 

19.30 

-13.1 

40.  00 

-14.2 
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The  initialization  procedure  is  to  take  the  temperature 
obtained  for  the  surface  in  the  previous  section,  and  use 
Table  4.2  to  arrive  at  a  guess  for  the  input  soil  profile. 

As  with  the  atmospheric  input  profile,  the  points  of  the 
input  soil  profile  need  not  be  co-aligned  with  the  soil 
grid. 


■ 
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CHAPTER  5 


RESULTS 


5-1  General 

In  this  chapter  the  results  from  the  computer 
simulations  of  Experiments  8,  9,  and  11  are  presented.  The 
results  are  discussed  in  relation  to  the  observed  data,  and 
observations  pertinent  to  the  model  are  also  included. 
Occasionally  recourse  is  made  to  observations  not  obtained 
at  the  Dawson  Bridge  site.  In  addition,  the  viability  of  a 
double  structure  inversion  in  the  valley  is  examined,  and 
the  differences  between  urban  and  rural  cooling  is  studied- 
Further,  sensitivity  to  certain  model  parameters  is 
considered.  In  each  experiment,  the  model  was  run  for  a 
period  of  12  hours  starting  at  local  noon.  The  grid 
employed  by  the  model  throughout  all  simulations  is  given  in 
Appendix  D.  Appendix  E  lists  the  parameters  which  were 
common  to  all  the  experiments,  and  Appendix  F  lists  the 
programme  for  the  model- 


' 
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5.2  The  fie suits  for  Exper iment  8 

5.  2.  1  General 

Experiment  8  was  conducted  on  27  June  1S78.  Conditions 
were  very  favourable.  No  significant  synoptic  or  mesoscale 
activity  existed  in  the  vicinity.  Scattered  fair-weather 
cumulus  prevailed  during  the  day,  and,  occasionally, 
scattered  altocumulus  occurred  during  the  evening.  The 
progression  of  wind  and  temperature  throughout  the  day  was 
almost  ideal.  At  noon  the  surface  winds  were  light.  Then, 
as  insolation  incited  buoyancy  in  the  boundary  layer, 
convective  overturning  mixed  the  stronger  winds  aloft 
downward.  A  mid-afternoon  maximum  of  6  ra/s  was  reached. 
Finally,  by  sunset  (approximately  1853  MST) ,  with  the 
decrease  of  surface  heating,  the  winds  again  became  light. 
Temperature  throughout  the  afternoon  was  essentially 
constant  although  small  variations,  likely  due  to  thermals, 
were  present.  A  maximum  of  26.5°C  (at  the  Municipal 
Airport)  was  reached  in  the  late  afternoon,  and  a  definite 
steady  temperature  decline  began  within  two  hours  of  sunset. 


5.2.2  The  Observations 

The  winds  and  temperatures  from  the  local  airports  are 
given  in  Table  5.1.  Valley  winds  are  given  Tables  5.2 
and  5.3.  Valley  temperatures  are  presented  in  Table  5.4. 
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Table  5.1  The  Temperatures  and  Winds  from  the  Local 
Airports  on  27  June  1978 


Airport 


Time 

(GMT) 

Municipal 

Namao 

International 

Temp 

(°C) 

Wind 
(°/  m/s) 

Temp 
(  °C ) 

W  ind 
(°/  m/s) 

Temp 

(°C) 

Kind 
(°/  m/s) 

1800 

22.  0 

300/3. 1 

21. 8 

290/6.2 

21.6 

290/4. 6 

1900 

22.  7 

28  0/2.  i 

2  2.  8 

280/6.2 

22.  4 

280/2.  6 

2000 

22.  9 

31 0/5. 7 

22.  7 

300/7. 7 

22.  6 

310/5. 7 

2100 

23.  1 

280/6. 2 

23.  7 

270/6. 7 

23.  3 

270/3. 1 

2200 

24.  1 

290/3. 1 

23.  S 

260/6.7 

24.  2 

2  70/4.  1 

2300 

24.  4 

270/2. b 

24.  2 

31 0/4. 6 

25.  1 

270/4. 1 

2400 

24.  8 

320/3.1 

24.  2 

310/2. 6 

25.2 

250/2. 6 

0000 

24.  6 

320/2. 6 

23.  9 

320/4.  6 

25.  4 

270/1.0 

otoo 

24.  0 

260/1. 5 

23.  9 

21 0/1.0 

24.2 

0 

0200 

23.  5 

270/1.0 

23.  6 

280/4. 0 

23.  1 

0  10/2.  1 

0300 

22.6 

0 

20.  8 

0 

18.  7 

0 

0400 

21.7 

0 

19.  4 

0 

16.  4 

0 

0500 

20.  2 

180/1. 0 

18.  8 

020/1.0 

14.  2 

120/1. 0 

0600 

16.3 

0 

17.  4 

0 

13.  6 

160/1.  0 

0700 

16.5 

0 

13.  0 

360/1. 5 

12.9 

150/1. 0 

0800 

15.  7 

0 

15.  4 

0 

12.  5 

310/3. 6 
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Table  5.  2 

The  Valley 

Winds 

(m/s)  on 

27  June 

1978 

Mean  Wind 
Station  11 

Mean  Wind 
Station  8 

Time 

(MDT) 

1.  70 

Height  (m) 

2.  84  5.84 

2.  64 1 

Height  ( m) 
2.  64 

1530-1600 

3.  95 

1600-1630 

4.  16 

1630-1700 

3.  30 

1 700-1730 

3.  70 

1730-1800 

1.  17 

1.  26 

1.  28 

1.21 

2.  96 

1800-1830 

0.  86 

0.  91 

0.  94 

0.  88 

3.  36 

1830-1900 

0.  89 

1.  00 

1.04 

0.  94 

2.  51 

1900-1930 

0.  76 

0.  89 

0.98 

0.  82 

2.  30 

1930-2000 

0.  52 

0.  55 

0.61 

0.  53 

2.  25 

2000-2030 

0.  61 

0.57 

0.  53 

0.59 

1.89 

2030-2100 

0.  43 

0.  35 

0.30 

0.  40 

1.  66 

2100-2130 

0.  41 

0.  42 

0.28 

0.41 

1.  25 

2130-2200 

0.  28 

0.  25 

0.  17 

0.27 

1.  05 

2200-2230 

0.20 

0.  32 

0.  29 

0.25 

0.  79 

2230-2300 

0.  11 

0.  33 

0.26 

0.  21 

0.  80 

2300-2330 

0.  15 

0.25 

0.23 

0.19 

1.  54 

2330-2400 

0.  09 

0.  08 

0.  07 

0.  09 

1.  26 

0000-0030 

0.  79 

^interpolated  linearly  from  1.70  m  and  2.84  m  winds. 


' 
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Table  5.3 

The  Slope  (Drainage) 
on  27  June  1978 

Wind  at 

Station  1 1 

lime 

Mean  Wind 

Time 

Mean  Wind 

(HOT) 

(m/s) 

(GMT) 

(m/s) 

1915 

-0.041 

2145 

0.  64 

20 

-0.041 

50 

0.  63 

25 

0.06 

55 

0.  64 

1930 

0.  06 

2200 

0.  66 

35 

0.38 

05 

0.  70 

40 

0.  28 

10 

0.  60 

1945 

0.1  7 

2215 

0.  56 

50 

0.  51 

20 

0.  66 

55 

0.38 

25 

0.  58 

2000 

0.15 

2230 

0.  71 

05 

0.59 

35 

Q.  51 

10 

0.21 

40 

0.  60 

2015 

0.68 

2245 

0.  70 

20 

0.42 

50 

0.  662 

25 

0.62 

55 

0.  622 

2030 

0.  77 

2300 

0.  572 

35 

0.  77 

05 

0.  53 

40 

0.49 

10 

0.  49 

2045 

0.71 

2315 

0.  43 

50 

0.  70 

20 

0.  62 

55 

0.54 

25 

0.  64 

2100 

0.85 

2330 

0.  45 

05 

0.72 

35 

0.  53 

10 

0.  66 

40 

0.  58 

2115 

0.  79 

2345 

0.  47 

20 

0.80 

50 

0.  36 

25 

0.71 

2130 

0.  86 

35 

0.81 

40 

0.  68 

iminus  sign  indicates  an  upslope  wind 
2linearly  interpolated  value. 


■ 
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Table  5.4  The  Calibrated  Valley  Temperatures  (°C)  on 


28  June  1978 

Station  8 

Station  11 

Time 

(MDT) 

Height 

(m) 

1.20  0.  36 

1.20 

Height 

(m) 

2.  60  5.  06 

7.  10 

9.  45 

1530  26.0 

45  25. 4 


1600  25.6 

15  25.9 

30  25.4 

45  25.6 

1700  26.6 


15 

25.3 

27. 

9 

25.  7 

25. 

6 

25.  4 

25. 

3 

30 

26.2 

26. 

9 

25.0 

24. 

5 

24.9 

25. 

1 

45 

26.  5 

27. 

1 

25.  6 

25. 

3 

25.  3 

25. 

1 

1800 

25.  8 

27. 

0 

25.6 

25. 

6 

25.5 

25. 

3 

15 

25.  8 

26. 

7 

25.  3 

25. 

3 

25.  7 

25. 

5 

30 

26.  0 

25. 

9 

24.6 

24. 

5 

25.3 

25. 
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5.2.3  The  Mode].  Parameters 

The  parameters  used  in  the  computer  simulation  have 
been  calculated  according  to  previously  detailed  procedures, 
and  are  listed  in  Table  5.5.  The  start-up  temperature 
profiles  for  the  atmosphere  and  soil  are  specified  in 
Tables  5.6  and  5.7.  The  2.64  meter  winds  required  by  the 
model  are  given  in  Table  5.8.  The  slope  friction  velocity 
was  calculated  according  to  equation  (3.4.4)  using  an  n  of 
0.27  determined  from  Table  5.2.  The  drainage  winds  required 
by  the  model  were  extracted  from  Table  5.3.  The  value  for  a 
missing  entry  was  linearly  interpolated  from  its  nearest 
neighbours.  Further,  in  the  interval  between  the 
termination  of  the  experimental  observations  and  the  model 
termination,  an  average  slope  wind  was  used. 


5.  2. 4  The  Model  Predictions 

In  this  section,  the  model  predictions  for  Experiment  8 
are  discussed.  Evidence  for  the  validity  of  the  calculated 
flux  density  profiles  is  presented;  the  characteristic s  of 
the  temperature  profiles  are  shown;  and  the  possibility  of  a 
double  structure  inversion  is  proffered. 

Generally,  the  magnitudes  and  trends  of  the  predicted 
flux  demsities  as  a  function  of  time  were  consistent  with 
the  literature  (see,  for  example,  Lettau  and  Davidson  1957)  » 
However,  proper  evaluation  of  the  flux  densities  produced  by 
the  numerical  model  was  difficult  since  it  was  not  possible 


. 
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Table  5.5  The  Model  Parameters  for  Experiment  8 


Paramete  r 

Progra  mine 
Name 

Val  ue 

Declination 

DEC 

23°  18.8* 

Equation  of  Time 

EOT 

2  min  58  sec 

Sun’s  Hour  Angle  at 

Model  Start-up 

KATO 

-1°  14.3* 

Maximum 

PHILIM 

1.069 

Slope  Wind  Start  Time 

TSLP 

18.58  hrs  M3T 

Surface  Pressure 

PB 

922  mb 

Average  Layer  Temperature 
throughout  Period 

TAV 

29  3.0  6  °C 

Atmospheric  Flux  Density 

FLWA 

2  65  W m~  2 

' 
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Table  5.6  The  Initial  Air  Temperature  Profile 
for  Experiment  8 


Height 

(m) 

Temperature 

(°K) 

2000. 

278.  96 

708. 

290.  16 

1. 20 

295.96 

1 . 00 

296.  16 

0.  46 

296.  76 

0.  22 

297.46 

0.  10 

298.16 

0.  05 

298. 86 

0.  00 

300. 52 

Table  5.7  The  Initial  Soil 
for  Experiment  8 

Temperature  Profile 

Depth 

Temperature 

(cm) 

(°K) 

0.  00 

300. 52 

0.  50 

298. 42 

1.  04 

296.32 

2.  15 

294. 12 

4.47 

291.82 

9.  28 

289. 62 

19.  30 

287. 42 

40.  00 

286.32 
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Table  5.8  The  Model  Hinds  (m/s)  for  Experiment  8 


Time 

Mean  Wind 

Mean  Wind 

Mean  Wind 

Friction 

Friction 

(GMT) 

(Municipal) 

(Observed ) 

(Filtered1 ) 

Velocity 

Velocity 

(Rim) 

( Slope) 

1800 

3.09 

1900 

2.06 

3.  60 

.258 

.159 

30 

4.  12 

.  295 

.182 

2000 

5.66 

4.63 

.332 

.  204 

30 

4.  99 

.358 

.  220 

2100 

6.  1  7 

5.34 

.383 

.236 

30 

4.95 

.  355 

.219 

2200 

4.  20 

4.55 

.326 

.201 

30 

4.  07 

.292 

.  180 

2300 

3.  28 

3.  56 

.257 

.158 

30 

3.35 

.  240 

.148 

0000 

3.  26 

3.  1 1 

.223 

.  137 

30 

2.98 

.214 

.  132 

0100 

2.  78 

2.85 

.204 

.126 

30 

2.  58 

.  185 

.114 

0200 

2.  52 

2.30 

.  165 

.  102 

30 

2.  01 

.  'i  44 

.  089 

0300 

1.61 

1. 72 

.123 

.076 

30 

1.42 

.  102 

.  06  3 

0400 

1.  03 

1.12 

.080 

.049 

30 

1.05 

.075 

.046 

0500 

0.  71 

0.97 

.070 

.043 

30 

0.97 

.070 

.  043 

0600 

1.  16 

0.96 

.069 

.042 

30 

1.01 

.072 

.  044 

0700 

1.002 

1. 05 

.0  75 

.046 

30 

1.03 

.074 

.045 

0800 

1.002 

1.00 

.072 

.  044 

30 

1.  00 

.072 

.  044 

0900 

1 . 002 

1.00 

.072 

.  04  4 

»a  three-point  Bartlett  filter  was  used. 

2calm  reading  assigned  a  value  of  1  m/s  as  per  section  3. 4 


■ 

63 


to  obtain  any  flux  measurements  in  any  of  the  river  valley 
experiments.  Nevertheless,  an  attempt  to  demonstrate  that 
the  model  flux  densities  are  reasonable  was  undertaken. 

Measured  net  radiation  and  solar  radiation  data9  for 
1977  are  available  from  Stony  Plain,  a  location 
approximately  35  kilometers  west  of  Edmonton.  Within 
certain  limitations,  the  assessment  of  the  predicted  solar 
flux  density  can  be  accomplished  using  this  data.  Errors 
arising  from  differences  in  location  were  presumed 
unimportant.  Errors  due  to  the  one  year  time  difference 
were  more  formidable.  Astronomically,  the  solar  radiation 
reaching  a  certain  point  at  a  particular  time  is  constant 
from  year  to  year.  Thus  any  difference  is  due  to 
atmospheric  effects  with  cloud  cover  in  particular  being  a 
major  factor.  Fortunately,  in  late  June  the  rate  of  change 
in  the  sun's  declination  is  a  minimum.  Thus,  it  can  be 
assumed  that  the  solar  radiation  is  similar  in  a  ten  day 
period  preceding  and  following  the  27  June  1977.  Further, 
it  is  assumed  that  the  variable  nature  of  cloud  cover  will 
ensure  that  each  of  the  24  hourly  observation  slots  will  be 
cloud-free  at  least  once  during  the  period  considered. 

Hence  a  good  estimation  of  the  solar  radiation  can  be  made 
by  taking,  for  each  hour,  the  maximum  value  recorded.  This 
procedure  was  carried  out,  and  the  results  plotted  in 


9This  data  is  published  in  the  Monthly  Radiation  Summary, 
ISSN  0027-0482,  by  the  Atmospheric  Environment  Service, 
Environment  Canada. 
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Figure  5.1.  The  solar  radiation  measured  at  Stony  Plain 
consisted  of  the  total  downward  direct  and  diffuse  solar 
flux;  no  measurement  of  reflected  solar  radiation  was  made. 
Thus  a  standard  reflection  of  20  per  cent  was  adopted.  The 
solar  radiation  calculated  by  the  model  is  also  plotted  in 
Figure  5.  1.  Given  the  magnitude  of  the  assumptions  made, 
the  agreement  is  quite  reasonable.  The  observed  maximum 
calculated  between  1100  and  1200  local  time  may  not  be  real. 
Significantly  lower  data  values  were  recorded  before  and 
after  this  maximum,  making  its  selection  suspect. 

A  similar  comparison  was  carried  out  for  the  net  flux 
of  radiation  reaching  the  earth's  surface.  The  results  are 
shown  in  Figure  5-2.  Since  the  net  flux  density  is  not 
necessarily  symmetric,  and  since  the  computer  simulation  was 
for  a  12  hour  period,  only  twelve  hours  of  the  Stony  Plain 
data  are  plotted.  The  results  are  not  out  of  line.  The  net 
flux  density  calculated  by  the  model  makes  no  allowance  for 
cloud  cover.  Thus  the  smaller  afternoon  amplitude  in  the 
calculated  observed  curve  is  likely  due  to  procedural 
inability  to  remove  completely  the  effect  of  cloud. 

Typically  summer  mornings  and  evenings  are  cloud-free,  while 
the  afternoons  develop  fair-weather  cumulus. 

The  model-predicted  latent  plus  sensible  heat  flux 
densities  to  the  atmosphere,  and  soil  heat  flux  density  to 
the  ground  are  given  in  Figures  5.3  and  5.4,  respectively- 
Unhappily,  the  Stony  Plain  data  does  not  include  these  flux 
measurements.  However,  the  verification  of  the  model  solar 
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Figure  5. 1  The  Assessment  of  the  Model-Predicted  Solar 
Flux  Density  for  Experiment  8 
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Figure  5.2  The  Assessment  of  the  Model-Predicted  Net  Flux 
Density  for  Experiment  8 
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Figure  5.3 


The  Model- Predicted  Heat  Flux  Density  to  the 
Atmosphere  for  Experiment  8 
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Figure  5.4  Ihe  Model-Predicted  Heat  Flux  Density  to  the 
Soil  for  Experiment  8 
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and  net  radiative  flux  densities,  and  the  general 
acceptability  of  the  magnitude  and  trend  in  the  remaining 
model- predicted  flux  densities  indicate  that  these  flux 
densities  are  sufficiently  reliable. 

The  evolution  of  the  model- predicted  temperature  fields 
for  both  the  rim  and  slope  between  1235  and  0035  MST  is 
presented  in  Figures  5.5  and  5.6.  In  order  to  more  closely 
scrutinize  the  results,  attention  will  be  focussed  on  two 
aspects  of  the  output:  the  variation  of  the  temperature  at 
a  particular  level  throughout  the  period;  and  the  variation 
of  the  temperature  in  the  vertical  at  a  particular  time. 

The  progression  of  predicted  temperature  with  time  for 
the  1  meter  level  over  the  slope  and  rim  is  shown  in 
Figure  5.7.  The  observed  screen  temperature  curves  are  also 
plotted  in  the  figure.  The  agreement  between  the  observed 
and  predicted  rim  temperatures  is  quite  good.  The  predicted 
curve,  although  generally  almost  one  degree  warmer,  is 
congruent  to  the  observed  trace,  and  the  times  of  maximum 
temperature  coincide.  Only  after  approximately  2300  I1DT 
does  the  predicted  trend  diverge  from  that  observed.  The 
slope  traces  do  not  agree  to  the  same  extent.  Initially  the 
agreement  is  good.  However,  once  the  slope  (drainage)  wind 
is  initiated,  the  trend  is  not  preserved,  and  by  2400  MDT, 
there  is  a  5  degree  difference  between  the  observed  and 
predicted  values.  The  inclusion  of  a  neutral  layer  in  the 
regime  of  the  slope  wind  tended  to  alleviate  this  problem. 

As  demonstrated  in  Figure  5.7,  the  trends  of  the  observed 
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figure  5.5  Ihe  Model- Predicted  Temperature  Field  at  the 
Bim  for  Experiment  8 
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Figure  5.6  The  Model- Predict ed  Temperature  Field  over 

the  Slope  with  the  Inclusion  of  the  Neutral  Flux 
Layer  for  Experiment  8 
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Figure  5a 7  Progression  of  the  1  meter  Temperature  with 

Time  for  Experiment  3 


and  predicted  curves  are  considerably  closer.  Although  the 
absolute  values  of  the  predicted  temperatures  are 
discouraging,  the  4.5  degree  difference  predicted  between 
the  rim  and  the  slope  stations  is  very  close  to  the  5  degree 
difference  observed.  Presumably  the  remaining  difference, 
and  the  actual  magnitude  of  the  temperatures,  could  be 
accounted  for  with  the  inclusion  in  the  model  of  both 
advection  and  a  more  sophisticated  mechanism  of  mixing  by 
the  slope  wind. 

The  vertical  temperature  profiles  are  considered  at 
four  times:  1615,  1800,  2000,  and  2200  MST.  The 
model-produced  profiles  for  the  slope  are  displayed  in 
Figure  5.8.  The  observed  profiles  below  10  metres  were 
actually  recorded.  The  'observed*  profiles  between  10  and 
90  meters  were  determined  from  the  Edmonton  City  Tower  data 
of  1977.  The  procedure  was  to  select  a  day  most  closely 
resembling  the  day  of  Experiment  8  and,  assuming  the 
temperature  gradients  between  10  and  90  meters  at  the 
appropriate  times  are  applicable,  use  them  to  calculate  an 
estimation  for  the  observed  profile.  Short  dashes  display 
these  crude  estimations  in  Figure  5.8. 

Comparison  of  the  computed  profiles  with  those  observed 
show  that  the  match  was  not  perfect.  However,  many 
essential  features  were  duplicated.  The  trends  in  the  1615 
MST  curves  were  in  good  agreement  even  though  the  computed 
curve  was  about  1  degree  too  warm.  Agreement  between  the 
1800  MST  profiles  was  poorer.  Again  the  computed  trace  is 
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figure  5.8  The  Vertical  Temperature  Profiles  over  the  Slope  for  Experimect 
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warmer  than  the  observed.  In  the  lowest  10  meters,  the 
observed  profile  indicated  that  the  formation  of  the 
inversion  was  well  under  way,  whereas  the  computed  profile 
was  just  approaching  neutrality.  Further,  the  observed 
curve  marked  a  much  more  unstable  trace  between  10  and 
50  meters  than  did  the  computed  profile-  However,  this  may 
have  been  a  result  of  the  extrapolation  used  to  arrive  at 
the  observed  curve  in  this  region.  Similarly,  comparison  of 
the  2000  and  2200  MS T  computed  and  observed  profiles  show 
that  the  trend  was  duplicated,  at  least  in  the  lower  levels, 
and  that  the  computed  curves  continued  to  be  warmer-  It  is 
disconcerting  that  the  trends  above  10  meters  did  not  align 
more  closely  with  observation.  This  is  difficult  to  account 
for;  and,  again,  the  methodology  of  the  extrapolation  is 
offered  as  the  explanation. 

The  distribution  of  heating  and  cooling  in  the  soil 
layer  follows  a  very  regular  course.  The  model  results  {see 
Figure  5.9)  are  exactly  those  expected  for  a  soil  layer  with 
a  constant  thermometric  conductivity  (see,  for  example. 
Sellers  1965).  During  periods  of  maximum  heating,  the 
uppermost  layers  in  the  soil  are  the  warmest.  The 
temperature  profile  is  essentially  linear  with  depth.  The 
lowest  layers  are  at  minimum  values.  Gradually,  the  surface 
heat  wave  propagates  downwards  through  the  soil.  Meanwhile, 
the  surface  layers  have  begun  to  cool.  Ey  late  evening  the 
result  is  a  zone  of  maximum  temperature  about  10  cm  deep 
with  cooler  soil  above  and  below.  The  temperatures  in  this 
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figure  5.9  The  Soil  Temperatures  for  Experiment  8 
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zone  are  not  nearly  as  high  as  the  mid-afternoon  maximums, 
for  the  temperature  wave,  as  it  moves  downwards,  dampens  by 
conduction  of  heat  to  the  cooling  surface  above  and  the 
cooler  depths  below  it.  Verification  of  the  model  soil 
predictions,  other  than  the  fact  that  proper  behaviour  of 
the  profiles  is  exhibited,  was  not  possible.. 

It  is  unrealistic  to  expect  that  the  profile  generated 
by  the  model  for  the  slope  would  be  legitimate  much  above 
rim  level.  Advection  and  mixing  between  the  more  extensive 
airmass  over  the  plain,  as  represented  by  the  rim  profile, 
and  the  air  over  the  valley  above  rim  level  should  ensure 
that  the  slope  profile  would  lose  its  identity.  If  the 
model  rim  and  slope  profiles  for  2200  MST  are  superimposed 
as  if  the  air  over  the  rim  had  drifted  over  the  valley  and 
mixed  only  in  a  10  meter  cushion  each  side  of  the  rim  level, 
then  a  double  structure  inversion  results  (Figure  5.10). 

This  structure  was  observed  by  Paterson  and  Hage  (1979),  and 
the  results  agree  qualitatively. 


5.3  Model  Sensitivity  to  Slope 

In  Experiment  8,  the  model-calculated  sunset  time  was 
18:53:29  MST.  This  time  was  determined  using  a  flat  uniform 
slope.  No  account  was  made  for  terrain  irregularities.  The 
predicted  temperatures  near  the  slope  surface  were 
consistently  warmer  than  observed.  A  possible  explanation 
for  this  trend  is  that  the  extinction  of  solar  radiation 
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Figure  5,10  The  Double-Structure  Inversion  for  Experiment  8 


actually  took  place  earlier  due  to  the  shading  from  trees. 
Alternatively,  it  could  be  that  the  slope  was  actually 
steeper  than  that  used  in  the  model.  To  test  the  viability 
of  these  explanations,  the  sensitivity  of  the  model  to  the 
value  used  for  the  inclination  of  the  slope  was  examined, 
run  for  Experiment  8  was  conducted  in  which  the  slope  was 
doubled.  A  comparison  of  the  1  meter  temperature 
progression  for  both  slopes  is  presented  in  Figure  5.11. 

The  slope  sunset  at  the  doubled  value  was  17:09:00  MST, 
approximately  an  hour  and  three-quarters  earlier.  The 
figure  shows  that  cooling  for  the  doubled  slope  preceded 
that  of  the  regular  slope  by  about  the  same  time.  The 
slightly  lower  initial  temperature  for  the  doubled  slope  is 
expected  in  light  of  the  decreased  insolation  received 
throughout  the  afternoon.  With  its  earlier  onset  of 
cooling,  the  steeper  slope  yields  a  temperature  more  than 
2.5  degrees  lower  by  midnight.  This  result  is  not 
unanticipated.  The  trends  of  the  two  curves  did  not  quite 
agree  around  1800  MDT.  This  is  due  to  the  fact  that  no 
allowance  was  made  to  start  the  model  slope  wind  at  an 
earlier  time.  Other  than  this  discrepancy,  the 
character ist ics  of  the  two  curves  matched.  Further,  as  can 
be  seen  from  Figures  5.7  and  5.11,  the  predicted  cooling 
rates  near  sunset  for  both  slopes  were  essentially  in 
agreement  with  those  observed.  Thus  it  is  concluded  that 
the  model  did  not  exhibit  any  undesired  sensitivity  to  the 
value  of  the  slope  used,  and  that  the  model  was  only 


' 
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Figure  5.11  The  Model  Sensitivity  to  Slope 
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sensitive  to  slope  in  that  it  determines  the  time  for  the 
onset  of  cooling.  Hence  the  higher  predicted  temperatures 
near  the  surface  can  indeed  be  due  to  an  improper 
consideration  of  the  physical  aspects  of  the  valley  and  to 
an  imprecise  slope  value.  It  is  also  surmised  that  the 
orientation  of  the  slope,  which  affects  the  onset  time  for 
cooling,  is  responsible  for  the  lower  temperatures  and 
stronger  inversions  recorded  over  the  east-  and  north-facing 
slopes. 


5.4  Model  Sensitivity  to  the  Soil  Dif f usivity 

To  determine  the  model  sensitivity  to  the  soil 
diffusivity,  the  run  for  Experiment  8  was  repeated  using  a 
thermal  conductivity  of  0.10  in  place  of  the 

previous  0.25  Wra-'K-1  value,  and  a  corresponding  thermal 
diffusivity  of  0.  6*10-7  m*2s~l  in  lieu  of  1.5*10~7 
The  resulting  1  meter  slope  temperature  curves  are  plotted 
in  Figure  5.12.  Figure  5-12  indicates  that  prior  to  sunset, 
the  model  experienced  very  little  sensitivity  to  the  soil. 
However,  once  the  atmosphere  became  stable,  and  the  sensible 
and  soil  heat  flux  densities  became  the  dominant  terms  in 
the  heat  balance  equation,  a  slight  sensivity  to  the  soil 
diffusivity  became  apparent..  By  2400  MDT,  close  to  0.  5  C° 
differences  in  the  1  meter  temperature  were  found  with  the 
40  per  cent  lower  soil  diffusivity. 

The  lower  temperatures  predicted  were  in  agreement  with 
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Figure  5.  12  The  Model  Sensitivity  to  the  Soil  Diffusivity 
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the  physics  of  the  model.  During  the  evening,  both  the 
sensible  and  soil  heat  flux  densities  are  directed  towards 
the  surface.  If  the  magnitude  of  the  soil  flux  density  is 
decreased,  then  the  heat  directed  to  the  surface  must  be 
decreased,  and  the  resulting  surface  temperature  lowered. 
By  conduction,  the  lower  surface  temperatues  yields  lower 
1  meter  temperatures. 

These  changes  in  evening  temperatures  indicate  some 
model  sensitivity  in  the  stable  regime  to  the  value  chosen 
for  the  soil  diffusivity. 


5.4  Model  Sensitivity  to  the  Friction  Velocity 

Figure  5.13  displays  the  1  meter  temperature  curves  for 
two  cases:  the  observed  friction  velocity  under  stable 
conditions  (about  .072  m/s),  and  one-quarter  of  this 
friction  velocity.  The  curves  are  divergent,  and  by 
midnight  more  than  an  1.5  degree  discrepancy  was  apparent. 

In  addition  to  predicting  a  decreased  cooling  rate  at  the 
1  meter  level  for  the  lower  friction  velocity,  the  model 
also  predicted  a  15  per  cent  more  intense  inversion  between 
the  1  meter  height  and  the  surface.  These  results  indicated 
that  the  model  is  quite  sensitive  to  wind,  and  that  a  major 
source  of  error  will  be  the  error  in  the  observed  winds. 

This  error  is  difficult  to  estimate  given  the  non-linear 
numerical  relationship  between  the  wind  and  the 


dif  fusivitie  s. 
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Figure  5.  13  The  Model  Sensitivity  to  the  Friction  Velocity 
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Although  the  sensitivity  to  wind  under  stable 
conditions  is  large,  it  is  not  unexpected.  In  the  stable 
case,  the  thermal  diffusivity  is  very  small,  and  it  is 
possible  that  not  all  the  heat  available  for  mixing  will  get 
mixed.  In  such  a  situation,  the  magnitude  of  the  thermal 
diffusivity  is  crucial.  Since  the  thermal  diffusivity  is 
essentially  proportional  to  the  friction  velocity,  the  model 
will  be  sensitive  to  this  wind  parameter.  Under  unstable 
conditions,  this  sensitivity  does  not  arise.  Presumably, 
the  mixing  is  always  sufficient  so  that,  regardless  of  a 
change  in  the  diffusivity,  all  of  the  available  heat  will  be 
completely  diffused. 

5.6  Comparison  between  Urban  and  Rural  Cooling 

In  Chapter  1,  it  was  stated  that  the  rural  temperatures 
and  nocturnal  cooling  rates  were  comparable  to  those  within 
the  valley.  Yet  the  mechanisms  involved  appear  quite 
distinct.  The  urban  environment  is  generally  windier,  and 
the  soil  diffusivity,  because  of  concrete,  is  higher. 
Further,  on  the  valley  slopes  there  are  the  drainage  winds 
and  earlier  extinction  of  solar  radiation.  To  ascertain 
whether  the  model  was  capable  of  accounting  for  the  rural 
temperature  decline,  an  evening  run  was  made.  The  run  used 
the  previous  values  of  wind  and  soil  diffusivity  for  the  rim 
station,  and,  at  the  rural  station,  used  one-guarter  of  the 
urban  wind,  and  one-fifth  of  the  soil  diffusivity. 
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Geiger  (1966)  indicates  that  the  diffusivity  for  soil  is 
about  an  order  of  magnitude  smaller  than  that  of  concrete. 
The  exact  contribution  of  concrete  to  the  effective  soil 
diffusivity  is  not  known,  so  the  rural  soil  difusivity  was 
arbitrarily  chosen  as  one-fifth  the  urban  value.  The 
results  for  the  1  meter  temperatures  are  displayed  in 
Figure  5.14.  It  is  encouraging  to  see  that  the  rural  curve 
was  indeed  lower  than  the  rim  trace.  However,  the 
difference  in  the  cooling  rate,  approximately  0.16  °C/hr, 
was  smaller  by  half  than  the  average  rate  of  0.32  °C/hr 
observed  for  the  same  period  (calculated  from  the 
International  and  Municipal  Airport  temperatures  in  Table 
5.1).  Further,  the  absolute  magnitudes  of  the  predicted 
temperatures  were  disappointing.  At  2400  MDT  the  predicted 
temperature  of  21.2  °C  was  much  larger  than  the  17.1  °C 
temperature  predicted  for  the  valley.  In  fact  it  was  only 
about  0.4  °C  lower  than  the  predicted  rim  temperature. 

Hence,  it  must  be  concluded  that,  to  explain  the  rural-urban 
difference  in  a  more  satisfactory  way,  the  model  needs 
additional  physics. 


Temperature  (C) 
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Time  (MST) 


Figure  5.  14  Comparison  between  Urban  and  Fural  Cooling 
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5.7  The  Results  for  Experiment  9 


5.7.1  General 

Experiment  9  was  conducted  on  04  July  1978,.  Conditions 
on  this  day  were  also  favourable  for  measurable  valley 
effects.  Initially  no  significant  synoptic  or  mesoscale 
activity  existed  in  the  area.  Scattered  fair-weather 
cumulus  prevailed  in  the  early  afternoon  and  developed  into 
scattered  towering  cumulus  by  late  afternoon.  In  the  early 
evening  all  convective  activity  halted,  and  the  sky 
condition  for  the  remainder  of  the  experiment  consisted  of  a 
few  tenths  of  stratocumulus  and  altocumulus,  with  some  thin 
cirrus.  During  the  afternoon  the  winds  were  initially 
light,  then  increased  to  a  mid-afternoon  peak  (3.4  m/s  at 
1600  MDT) ,  after  which  they  gradually  decreased,  becoming 
light  just  before  sunset.  However,  at  about  2200  MDT  the 
wind  began  to  increase  again  reaching  an  overnight  maximum 
of  3.6  m/s.  Further,  after  the  experiment  had  terminated, 
the  stratocumulus  thickened  up  to  a  broken  condition,  and  at 
midnight  the  Municipal  Airport  recorded  a  light  shower.  The 
exact  nature  of  this  mesoscale  activity  is  not  known.  Thus 
conditions  for  Experiment  9  were  not  as  ideal  as  for 
Experiment  8  for  three  reasons:  firstly,  the  stronger 
afternoon  convective  activity  indicated  an  afternoon  latent 
heat  flux  density  greater  than  the  ideal  case  assumed  by  the 
model;  secondly,  the  mid-evening  increase  in  wind  resulted 


. 


in  a  larger  diffusivity;  and,  thirdly,  the  increased 
mid-evening  cloud  cover  resulted  in  a  greater  than  assumed 
atmospheric  long-wave  flux  density- 
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5-7.2  The  Observations 

The  winds  and  temperatures  from  the  local  airports  are 
given  in  Table  5.9.  Valley  winds  are  displayed  in 
Tables  5-10  and  5.11.  Table  5-12  presents  the  valley 
temperatures.  In  addition  to  the  complicating 
meteorological  factors  mentioned  in  the  previous  section, 
there  was  an  observational  difficulty.  The  rim  station 
anemometer  was  periodically  unreliable  during  the  evening 
with  a  mechanical  problem.  Thus  the  valley  data  had  to  be 
augmented  further  with  airport  data. 


5.7.3  The  Model  Parameters 

The  computer  simulation  for  Experiment  9  used  a  set  of 
parameters  similar  to  Experiment  8.  These  parameters  were 
calculated  in  an  identical  fashion,  and  are  given  in 
Table  5.13.  The  start-up  profiles  for  the  atmosphere  and 
soil  are  specified  in  Tables  5.14  and  5.15.  The  2.64  meter 
winds  required  by  the  model  are  given  in  Table  5.16.  An  n 
of  0.42  was  determined  from  Table  5.10,  and  Table  5.11  also 
lists  the  drainage  winds  needed  by  the  model..  For  missing 
values,  interpolation  and  extrapolation  procedures  imitated 
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Table  5.9  The  Temperatures  and  Hinds  from  the  Local 
Airports  on  04  July  1978 


Airport 


Municipal  Namao  International 


Time 

(GMT) 

Temp 

(°C) 

Wind 
(°/  m/s) 

Temp 

(°C) 

Wind 
(°/  m/s) 

Temp 

(°C) 

Wind 
(°/  m/s) 

1800 

24.3 

290/1. 5 

23.  6 

270/3. 1 

24.  2 

270/2. 1 

1900 

25.  7 

280/2. 1 

25.  2 

270/2. 1 

25.  9 

230/1.  5 

2000 

25.5 

270/3.1 

26.  3 

260/1.5 

26.6 

250/3.6 

2100 

26.  6 

220/2. 6 

24.  2 

270/4. 1 

27.  9 

210/2. 6 

2200 

27.2 

360/1. 0 

25.  8 

270/2.6 

27.  2 

220/1. 0 

2300 

26.9 

260/2. 6 

2  7.  7 

360/4. 6 

26.9 

300/3.  1 

0000 

27.1 

300/2. 1 

27.  6 

310/1.5 

27.  2 

290/2.6 

0100 

26.  b 

340/1. 5 

2  7.  0 

300/2.  6 

26.  8 

280/1. 5 

0200 

26.3 

290/2.  1 

26.  1 

310/1.5 

25.8 

320/2.  1 

0300 

24.  1 

330/1. 0 

24.  9 

0 

25.  2 

330/1. 0 

0400 

20.  9 

330/1.0 

20.  1 

0 

20.  8 

250/1.0 

0500 

20.  2 

070/1.5 

18.  6 

030/2.  1 

19.  1 

250/1. 0 

0600 

19.6 

050/4. 1 

1  7.  4 

020/4.1 

19.3 

0  10/5.  7 

0700 

17.0 

020/2. to 

16.  6 

030/5.  1 

18.  1 

020/4.  1 

0800 

16.  3 

070/3.6 

15.  4 

070/2. 6 

16.8 

030/3. 1 

0900 

15.  5 

030/2. 6 

14.  0 

020/1.  5 

16.  1 

030/2. 6 

- 
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Table  5.  10 

The  Valley 

Hinds 

(m/s)  on 

04  July 

1978 

Mean  Hind 
Station  11 

Mean  Wind 
Station  8 

Time 

(MDT) 

1.  70 

Height  (m) 

2.  84  5.  84 

2.  64i 

Height  (m) 

2.  64 

1500-1530 

2.  48 

1530-1600 

2.  39 

1600-1630 

1.  52 

1630-1700 

1.  1 1 

1.  17 

1.20 

1.14 

1.  79 

1700-1 730 

0.  96 

1. 09 

1.  14 

1.02 

1. 48 

1730-1800 

1.  20 

1.27 

1. 29 

4.  23 

2.  27 

1800-1830 

0.  93 

1.02 

1.36 

0.97 

1.  40 

1830-1900 

0.  80 

0.  86 

0.  S3 

0.  83 

1.  51 

1900-1930 

0.  62 

0.  72 

0.78 

0.66 

0.  71 

1930-2000 

0.  57 

0.  58 

0.55 

0.  57 

1.  07 

2000-2030 

0.  54 

0.  50 

0.46 

0.52 

0.  44 

2030-2100 

0.  53 

0.  38 

0.  43 

0.46 

0.  052 

2100-2130 

0.  33 

0.  31 

0.  24 

0.  32 

0.  022 

2130-2200 

0.  23 

0.  19 

0.09 

0.21 

0.  28 

2 200-2230 

0.  43 

0.  50 

0.  58 

0.46 

2.  34 

2230-2300 

0.  41 

interpolated  linearly  from  1,70  m  and  2.84  m  winds 
2suspect;  not  used. 
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Table  5.11 

The  Slope  (Drainage) 
on  04  July  1978 

Wind  at 

Station  11 

Time 

Mean  Wind 

Time 

Mean  Wind 

(MDT) 

(m/s) 

(GMT) 

(m/s) 

1835 

-0.291 

2035 

0.  65 

40 

-0. 56i 

40 

0.62 

1845 

0.02 

2045 

0.  58 

50 

-0.411 

50 

0.  63 

55 

-o. i b i 

55 

0.  71 

1900 

0.06 

2100 

0.  68 

05 

0.  02 

05 

0.  74 

10 

0.19 

10 

0.  68 

19  15 

-0. 04i 

21  i  5 

0.  68 

20 

0.22 

20 

0.  72 

25 

0.38 

25 

0.  59 

1930 

0.00 

2130 

0.  57 

35 

0.  42 

35 

0.  62 

40 

0.62 

40 

0.  73 

1945 

0.  79 

2145 

0.  54 

50 

0.26 

50 

0.  54 

55 

0.  34 

55 

0.  66 

2000 

0.  15 

2200 

0.  49 

05 

0.38 

05 

0.  30 

10 

0.47 

10 

0.57 

2015 

0.  32 

2215 

-o.  in 

20 

0.56 

20 

0.  19 

25 

0.  53 

25 

0.  45 

2030 

0.56 

iminus  sign  indicates  an  upslope  wind 
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Table  5.  12  The  Calibrated  Valley  Temperatures  (°C)  on 
04  July  1978 


Time 
( M  D  T ) 

Station  8 

Station  11 

Height 

(m) 

1. 20 

1.  20 

2.  44 

Height 

(m) 

4.  82 

7.  10 

9.45 

1600 

29.  3 

15 

27.9 

30 

29.0 

29.  4 

27.  7 

27.4 

28.  1 

27.  6 

45 

28.7 

29.  3 

28.  2 

28.  0 

28.3 

28.4 

1700 

29.  1 

29.  5 

28.  7 

28.  0 

28.  2 

27.  8 

15 

29.6 

29. 4 

28.2 

28.  0 

28.4 

28.  3 

30 

29.2 

29.  4 

28.  2 

28.  1 

28.  6 

28.  2 

45 

29.4 

29.3 

28.  5 

28.2 

28.4 

28.  2 

1800 

29.3 

29.  0 

28.  6 

28.  9 

29.  0 

28.5 

15 

29.2 

29.  0 

28.  0 

28.  1 

28.  3 

28.  0 

30 

28.9 

28.  7 

27.  6 

28.  0 

28.  1 

28.  3 

45 

29.1 

28.  1 

28.  2 

28.3 

28.  6 

28.2 

1900 

29.  2 

27.  4 

28.  4 

28.  2 

28.  2 

28.  1 

15 

29.5 

26.  4 

27.4 

27.4 

27.  8 

27.  7 

30 

28.7 

25.  6 

26.  0 

26.7 

27.  7 

27.  1 

45 

28.0 

24.0 

26.2 

26.  7 

27.  0 

27.2 

2000 

28.  4 

24.  7 

23.  5 

25.  8 

27.  1 

27.  2 

15 

28.  0 

23.  2 

24. 9 

26.6 

26.6 

26.  8 

30 

28.2 

22.  5 

22.  2 

23.4 

23.  9 

25.  2 

45 

27.4 

21. 9 

23.  6 

25.8 

25.6 

26. 3 

2100 

26.  5 

22.  0 

22.3 

25.  7 

25.  6 

26.  1 

15 

25.8 

21. 0 

21.5 

24.3 

25.  4 

25.  9 

30 

25.4 

20.  1 

21.  2 

22.6 

24.  1 

24.  9 

45 

25.3 

19. 3 

22.  0 

23.6 

23.4 

24.  7 

2200 

24.  3 

19.  6 

22.2 

22.  8 

22.  8 

23.  7 

15 

22.  2 

20.  8 

22.  0 

22.3 

22.  3 

22.  7 

30 

21.5 

20.  5 

21.5 

21.3 

21.3 

21.  6 

45 

21.3 

2300 

21.  1 

■ 

Table  5.13  The  Model  Parameters  for  Experiment  9 


Paramete  r 


Programme 

Name 


Value 


DEC 


22°  51.3 


Declination 

Eguation  of  Time 

Sun's  Hour  Angle  at 
Model  Start-up 

Maximum  ^ 

Slope  Wind  Start  Time 

Surface  Pressure 

Average  Layer  Temperature 
throughout  Period 

Atmospheric  Flux  Density 


EOT 

4  min  19  sec 

HATO 

-1°  54.7* 

PHIL  III 

0.740 

TSLP 

18.58  hrs  MST 

PB 

922  mb 

TAV 

295-  34  °C 

FLWA 

260  Wm-2 
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Table  5.14  The  Initial  Air  Temperature  Profile 
for  Experiment  9 


Height 

(m) 

Temperature 

<°K) 

2000. 

279. 83 

702. 

292.  16 

1. 20 

298. 76 

1.00 

298.96 

0.  46 

299. 56 

0.  22 

300.26 

0.  10 

300. 96 

0.  05 

301. 66 

0.  00 

303. 32 

Table  5.15  The 
for 

Initial  Soil 
Experiment  9 

Temperature  Profile 

Depth 

Temperature 

(cm) 

(°K) 

0.  00 

303.32 

0.  50 

301.  22 

1.  04 

299.12 

2.  15 

296. 92 

4.47 

294.62 

9.  28 

292.  42 

19.  30 

290. 22 

40.  00 

289.12 

96 


Table  5. 16  The  Model  Winds  (m/s)  for  Experiment  9 


Time 

(GMT) 

Mean  Wind 
(Municipal) 

Mean  Wind 
(Observed) 

Mean  Wind 
(Filtered1) 

Friction 

Velocity 

(Him) 

Friction 

Velocity 

(Slope) 

1800 

1.54 

1900 

2.06 

2.  23 

.160 

.153 

30 

2.40 

.  172 

.165 

2000 

3.09 

2.57 

.  184 

.  1  76 

30 

2.  79 

.200 

.'191 

2100 

2.57 

3.  00 

.215 

.206 

30 

2.  78 

.  200 

.191 

2200 

3.  35 

2.56 

.  184 

.  176 

30 

2.  54 

.  183 

.  175 

2300 

1.  77 

2.52 

.181 

.  1  73 

30 

2.  17 

.156 

.149 

0000 

2.  43 

1.81 

.  130 

.  124 

30 

1.81 

.  130 

.  124 

0100 

1. 22 

1.80 

.129 

.124 

30 

1.57 

.  112 

.108 

0200 

1.  74 

1.33 

.095 

.091 

30 

1.30 

.  093 

.  089 

0300 

1.  03 

0.  012 

1.27 

.091 

.087 

30 

1.24 

.  089 

.  085 

0400 

1.03 

0.  542 

1.20 

.  086 

.  082 

30 

1.  72 

.  123 

.118 

0500 

1 . 54 

0.  052 

2.23 

.  160 

.153 

30 

2.49 

.  178 

.171 

0600 

4.  12 

2.  74 

.  196 

.188 

30 

3.  09 

.  221 

.212 

0  700 

2.57 

3.  43 

.246 

.236 

30 

3.  17 

.  228 

.218 

0800 

3.  60 

2.91 

.209 

.200 

30 

2.  83 

.203 

.  194 

0900 

2.57 

2.  74 

.196 

.188 

ia  three-point  Bartlett  filter  was  used. 
2suspect;  not  used. 
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those  of  Experiment  8, 


5.7.4  The  Model  Predictions 

This  section  discusses  the  model  predictions  for 
Experiment  9.  The  main  intention  is  to  demonstrate  that  the 
computer  model  was  not  tuned  specifically  to  Experiment  8. 
The  forecast  temperature  fields  between  1235  and  0035  MST 
for  the  rim  and  the  slope  are  shown  in  Figures  5.15  and 
5.16.  It  is  seen  that  the  model  behaved  properly:  the 
general  features  were  consistent  with  the  observations.  The 
progression  of  the  1  meter  temperatures  with  time  is  given 
Figure  5.17.  Again  there  is  a  reasonable  correlation 
between  the  predicted  and  observed  trends.  However,  the 
absolute  values  predicted  were  discouraging.  Near  2200  MDT 
the  observed  temperature  increased.  Increased  mixing, 
attributed  to  the  increase  in  observed  wind  at  this  time,  is 
believed  responsible  for  the  increase  in  the  observed 
temperature.  Greater  mixing  transports  the  warmer  air  in 
the  inversion  down  to  the  surface.  Because  of  the  effect  of 
the  three-point  smoothing  employed,  the  model  winds  were  not 
increased  until  after  2300  MDT.  Thus  the  predicted 
temperature  increase  by  the  model  was  delayed. 


" 
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Figure  5.15  The  Model-Predicted  Temperature  Field  at  the 

Eim  for  Experiment  9 
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Figure  5.16  The  Model-Predicted  Temperature  Field  over 

the  Slope  with  the  Inclusion  of  the  Neutral 
Surface  Layer  for  Experiment  9 
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Figure  5*17  Progression  of  the  1  meter  Temperature  with 
time  for  Experiment  9 
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5.8  The  Results  for  Experiment  JlJ 

5.8.1  General 

Experiment  11  was  conducted  on  28  August  1978.  Once 
again,  weather  conditions  on  the  day  were  initially 
favourable.  The  early  afternoon  saw  scattered  cumulus 
clouds  with  fairly  brisk  winds  hovering  near  the  5  m/s  mark. 
By  mid-afternoon  cumulonimbus  clouds  had  developed  in  the 
area  and  by  1700  MDT  there  was  a  line  of  thundershowers  to 
the  north.  At  approximately  1815  MDT  the  wind,  which  had 
been  slowly  decreasing,  shifted  from  the  northwest  to  the 
northeast  and  picked  up  to  4.7  m/s  in  strength. 
Simultaneously,  the  temperature  dropped  nearly  2  degrees  in 
about  10  minutes.  Virga  was  observed  in  the  distance. 
However,  by  1900  MDT,  this  mesoscale  disturbance  appeared  to 
have  moved  completely  through  the  valley,  and  the  usual  rate 
of  decrease  in  both  wind  speed  and  temperature  resumed. 
However,  at  2300  MST  the  wind  again  picked  up  before 
dropping  to  near  calm  values  for  the  remainder  of  the 
period.  Thus  conditions  for  Experiment  11  were  not 
completely  ideal,  again  for  three  reasons:  firstly,  the 
significant  convective  activity  indicated  an  afternoon 
latent  heat  flux  density  greater  than  the  ideal  case  assumed 
by  the  model;  secondly,  the  mesoscale  frontal  passage 
introduced  a  2  C°  drop  which  the  model  can  not  take  into 
account;  and,  thirdly,  a  late  evening  increase  in  the  wind 
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resulted  in  a  lar ger-than-ideal  diffusivity  for  that  period. 


5.8.2  The  Observations 

The  winds  and  temperatures  from  the  local  airports  are 
given  in  Table  5.  17.  Valley  winds  are  displayed  in 
Tables  5.18  and  5.19.  Table  5.20  presents  the  valley 
tem  pera tures. 


5.8.3  The  Model  Parameters 

A  set  of  parameters  similar  to  those  of  the  previous 
two  experiments  was  used.  These  are  listed  in  Table  5.21. 
The  start-up  profiles  for  the  air  and  soil  are  given  in 
Tables  5.22  and  5.23.  The  2.64  meter  winds  required  are 
given  in  Table  5.24.  An  n  of  0.25  was  determined  from 
Table  5.18.  Table  5.19  specifies  the  drainage  (slope)  winds 
used  by  the  model.  For  missing  values,  the  previous 
interpolation  and  extrapolation  procedures  were  used- 


5.8.4  The  Model  Predictions 

In  this  section  the  model  predictions  for 
Experiment  11,  like  Experiment  9,  are  very  briefly 
presented.  Again  the  main  intention  is  to  show  the  model 
capablity  under  different  situations.  The  forecast 
temperature  fields  between  1235  and  0035  MST  for  the  rim  and 
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Table  5,17  The  Temperatures  and  Winds  from  the  Local 
Airports  on  28  August  1978 


Airport 

Municipal 

Namao 

International 

Time 

(GMT) 

Temp 

<°C) 

Wind 
(°/  m/s) 

Temp 

(°C) 

Wind 
(°/  m/s) 

Temp 

(°C) 

Wind 
(°/  m/s) 

1800 

17.9 

310/4,.  6 

18.4 

290/6.2 

17.5 

280/2. 6 

1900 

19.4 

280/5. 1 

19.  1 

290/5. 7 

G8.  7 

19.  1 

280/6. 6 

2000 

19.9 

290/5. 1 

1  8.  2 

290/5. 7 
G9.3 

19.  7 

290/4. 1 

2100 

19.  7 

300/4. 6 

19.  6 

290/6.2 
G10. 3 

19.6 

300/5. 1 

2200 

20.  4 

290/5.  7 

20.  5 

290/5. 7 
G8.7 

20.  2 

290/5. 7 

2300 

2  0.  4 

330/4. 1 

20.  6 

300/4. 1 

20.2 

320/6.  2 

0000 

20.  5 

31 0/4.  1 

20.  4 

350/5. 7 

G 1 1.  3 

20.4 

290/4. 1 

0100 

18.  1 

020/3. 6 

16.  9 

330/3. 1 

19. 4 

320/3. 6 

0200 

17.0 

010/2.  1 

16.  4 

310/1. 5 

16.7 

040/1.  5 

0300 

15.  2 

01 0/1.0 

14.  2 

280/1.0 

13.3 

0 

0400 

15.  2 

040/3.  1 

12.  1 

360/4.  1 

10.  5 

0 

0500 

12.  1 

070/1. 5 

11.8 

050/1.5 

9.9 

0 

0600 

11.  8 

0 

10.  6 

050/1. 0 

8.  1 

0 

0700 

11.6 

0 

11.  1 

0 

7.  4 

190/1. 0 

0800 

10.  3 

0 

9.  9 

0 

7.  0 

190/1. 5 

0900 

10.3 

0 

11.3 

0 

6.9 

200/2. 1 

Table  5.18  The  Valley  Winds  (m/s)  on  28  August  1978 


Time 

(MDT) 

Mean  Wind 
Station  11 

Mean  Wind 
Station  8 

1.  70 

Height 
2.  84 

(m) 

5.84 

2.  64 1 

Height  (m) 
2.  64 

1530-1600 

4.  47 

1600-1630 

4.  80 

1630-1700 

4.  32 

1 700-1 730 

4.  22 

1730-1800 

1 . 20 

1. 33 

1.37 

1. 26 

3.33 

1800-1830 

1.  02 

1.  08 

1.  15 

7.0  5 

2.  58 

1830-1900 

1.  14 

1.33 

1.54 

1 .22 

3.  86 

1900-1930 

1.  07 

1.  38 

1.  73 

1.21 

4.  72 

1930-2000 

0.  70 

0.  93 

1.  19 

0.  80 

2.  98 

2000-2030 

0.  46 

0.  43 

0.50 

0.  45 

1. 86 

2030-2100 

0.  57 

0.47 

0.37 

0.53 

2.  04 

2100-2130 

0.  43 

0.  41 

0.  44 

0.42 

1.  77 

2130-2200 

0.23 

0.  35 

0.45 

0.28 

0.  98 

2200-2230 

0.  31 

0.  36 

0.37 

0.  33 

1.  62 

2230-2300 

0.  36 

0.  39 

0.37 

0.37 

2.  55 

2300-2330 

0.  51 

0.  48 

0.58 

0.50 

2.  80 

2330-2400 

0.  35 

0.  42 

0.45 

0.38 

1. 51 

0000-0030 

0.  15 

0.  08 

0.05 

0.  12 

1.  92 

0030-0100 

1.  02 

linterpolated  linearly  from  1.70  m  and  2.84  m  winds 


• 
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Table  5.19  The  Slope  (Drainage)  Kind  at  Station  11 
on  28  August  1978 


Time 

(MDT) 

Mean  Wind 
(m/s) 

Time 

(GMT) 

Mean  Kind 
(m/s) 

1830 

0.00 

2115 

0.  17 

35 

-0. 071 

20 

0.  12 

40 

0.00 

25 

0.  1  7 

1845 

0.  00 

2130 

0.  15 

50 

0.12 

35 

0.  21 

55 

0.  00 

40 

0.  19 

1900 

0.04 

2145 

0.  19 

05 

-0.041 

50 

0.  20 

10 

0.12 

55 

0.  05 

1915 

0.  12 

2200 

0.  06 

20 

0.02 

05 

0.  30 

25 

0.  12 

10 

0.  23 

1930 

0.21 

2215 

0.31 

35 

0.  33 

20 

0.  16 

40 

0.22 

25 

0.  10 

1945 

0.  32 

2230 

0.  23 

50 

0.45 

35 

0.  30 

55 

0.  49 

40 

0.  27 

2000 

0.41 

2245 

0.41 

05 

0. 52 

50 

0.  40 

10 

0.48 

55 

0.  42 

2015 

0.  59 

2300 

o 

. 

00 

20 

0.51 

05 

0.  44 

25 

0.  48 

10 

0.  44 

2030 

0.44 

2315 

0.  36 

35 

0.  43 

20 

0.  33 

40 

0.31 

25 

0.  42 

2045 

0.  29 

2330 

0.  33 

50 

0.33 

35 

0.  28 

55 

0.36 

40 

0.  23 

2100 

0.  32 

2345 

0.  23 

05 

0.  25 

50 

0.  32 

10 

0.22 

55 

0.  26 

iniinus  sign  indicates  an  upslope  wind 
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Table  5. 20  The  Calibrated  Valley  Temperatures  (°C)  on 
27  August  1978 


Station  8 

Station  1 1 

Time 

(MOT) 

Height 

(m) 

1.  20 

0.  12 

1.20 

Height 

(m) 

3.  72 

6.  89 

14.  73 

1515 

19.8 

30 

20.  5 

45 

20.  7 

1600 

21.0 

15 

20.9 

30 

21.  1 

45 

21.  1 

1700 

21.2 

22.  3 

23.  3 

24.  4 

23.  5 

15 

21 . 1 

20.  3 

22.  1 

20.  5 

20.4 

18.8 

30 

21.  1 

19.  7 

21.6 

20.  6 

20.  6 

20.  7 

45 

21.6 

19.  8 

21 . 1 

20.6 

20.  7 

21.  2 

1800 

21.4 

18.  7 

20.  8 

20.  7 

21.2 

21.  3 

15 

21.4 

1  7.  6 

20.0 

20. 5 

21 . 0 

21.3 

30 

19.  3 

18.  2 

18.  6 

18.  7 

18.  7 

19.  0 

45 

18.  5 

18.  2 

18.0 

18.  7 

16.  6 

18.  9 

1900 

18.6 

7  8.  4 

18.  0 

18.  9 

18.  9 

19.  2 

15 

18.6 

17.9 

1  7.9 

18.  6 

18.8 

18.  8 

30 

17.7 

17.  2 

17.2 

18.  3 

18.  7 

19.  0 

45 

17.  9 

15.  1 

16.3 

17.  8 

18.4 

18.  7 

2000 

17.9 

13.  0 

15.  0 

17.  0 

77.  1 

18.  6 

15 

17.2 

12.  6 

14.  1 

16.8 

18.3 

19.1 

30 

16.  8 

13.  3 

13.9 

7  5.  6 

18.  1 

18.  7 

45 

16.  8 

12.  2 

13.7 

15.9 

17.0 

17.  5 

2100 

16.4 

11.  7 

13.  0 

14.  6 

75.  6 

16.  6 

15 

16.2 

11.6 

12.  6 

14.2 

15.5 

16.  9 

30 

16.  0 

10.  6 

12.  4 

14.  0 

14.  7 

16.  4 

45 

15.  8 

10.  8 

12.  0 

13.9 

14. 8 

15.  9 

2200 

15.  8 

misg 

12.  3 

misg 

misg 

misg 

15 

15. 1 

10.  4 

11.8 

13.  4 

14.6 

15.3 

30 

1  4.  4 

9.  0 

11.7 

7  2.  8 

13.  6 

14.  6 

45 

13.9 

10.  2 

11.3 

13.  0 

14.0 

14.  6 

2300 

13.  5 

8.  7 

11.2 

13.  7 

14.  4 

14.  5 

15 

13.2 

9.  1 

10.  6 

12.5 

13.5 

14.2 

30 

13.  0 

8.  7 

10.2 

12.  8 

13.  4 

14.  0 

45 

12.  6 

8.  3 

10.1 

11.9 

12.  5 

13.  5 

0000 

12.5 

8.  1 

10.  u 

11.5 

7  1.7 

12.  7 

15 

12.4 

- 
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Table  5,21  The  Model  Parameters  for  Experiment  11 


Parameter 

Programme 

Name 

Value 

Declination 

DEC 

9°  37.7’ 

Equation  of  Time 

EOT 

1  min  11  sec 

Sun’s  Hour  Angle  at 

Model  Start-up 

HATO 

1 

o 

o 

—a 

• 

m 

Maximum  <f>h 

PHILIM 

1.824 

Slope  Wind  Start  Time 

TSLP 

18.58  hrs  MST 

Surface  Pressure 

PB 

930  mb 

Average  Layer  Temperature 
throughout  Period 

TAV 

286.96  °C 

Atmospheric  Flux  Density 

FLWA 

245  Wm-2 

’ 


■ 
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Table  5.22  The 
for 

Initial  Air  Temperature  Profile 

Experiment  11 

Height 

(m) 

Temperature 

<°K) 

2000. 

275. 49 

1497. 

278.  96 

983. 

282. 76 

708. 

284.96 

1.20 

292. 85 

1.  00 

293. 05 

0.  46 

293.65 

0.  22 

294.35 

0.  10 

295. 05 

0.  05 

295.  75 

0.  00 

297. 41 

Table  5.23  The 
for 

Initial  Soil  Temperature  Profile 

Experiment  11 

Depth 

Temperature 

(cm) 

(°K) 

0.  00 

297.  41 

0.  50 

295.31 

1. 04 

293.  2*1 

2.  15 

291.01 

4.  47 

288.  71 

9.  28 

286.51 

19.  30 

284.31 

40.  00 

283. 21 
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Table 

5.  24  The  Model  Winds 

(m/s)  for  Experiment  11 

Time 

Mean  Wind  Mean  Wind 

Mean  Wind 

Friction 

Friction 

(GMT) 

(Municipal)  (Observed) 

(Filtered1 2) 

Velocity 

Velocity 

(Bim) 

( Slope) 

1800 

4.63 

1900 

5.  14 

5.  06 

.363 

.  207 

30 

5.  14 

.368 

.210 

2000 

5.  14 

5.06 

.  363 

.207 

30 

4.89 

.351 

.200 

2100 

4.63 

4.  66 

.334 

.190 

30 

4.  63 

.332 

.189 

2200 

4.  80 

4.  53 

.  325 

.185 

30 

4.45 

.319 

.182 

2300 

4.  22 

3.  96 

.  284 

.  162 

30 

3.  38 

.242 

.138 

0000 

2.  58 

3.  26 

.  234 

.133 

30 

3.  72 

.267 

.  152 

0100 

4.  72 

3.  85 

.276 

.  157 

30 

3.  19 

.229 

.131 

0200 

1.  86 

2.29 

.  1  64 

.  093 

30 

1.89 

.  135 

.077 

0300 

1.  77 

1.60 

.  115 

.066 

30 

1 . 46 

.105 

.060 

0400 

1.  62 

1.72 

.  123 

.070 

30 

2.32 

.  166 

.095 

0500 

2.  80 

2.  29 

.164 

.  093 

30 

2.  08 

.  149 

.085 

0600 

1.  92 

1.48 

.  106 

.060 

30 

1.31 

.094 

.054 

0700 

1.002 

1.01 

.072 

.  041 

30 

1.00 

.072 

.  041 

0800 

1. 002 

1.00 

.072 

.  041 

30 

1.00 

.072 

.  041 

0900 

1.  002 

1.  00 

.072 

.  041 

la  three-point  Bartlett  filter  was  used. 

2calm  reading  assigned  a  value  of  1  m/s  as  per  section  3. 4. 
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the  slope  are  shown  in  Figures  5-18  and  5,19.  The  results 
again  demonstrate  that  the  main  features  were  consistent 
with  the  observations-  For  closer  inspection,  the  evolution 
of  the  1  meter  temperatures  is  given  in  Figure  5.20.  It  is 
seen  that  when  the  2  degree  cooling,  which  occurred  with  the 
mesoscale  cold  front,  is  removed  from  the  observations,  the 
relationship  between  the  predicted  and  observed  curves  was 
very  close  to  those  of  Experiments  8  and  9-  After  2300  MDT, 
the  predicted  rim  curve  depicts  a  decrease  in  the  cooling 
rate,  and  the  predicted  slope  curve  shows  a  temperature 
increase.  This  is  consistent  with  the  increased  model  winds 
used  for  this  period.  Why  this  warming  was  not  more 
strongly  reflected  in  the  observations  (a  decrease  in  the 
cooling  rate  only  is  detected)  is  not  known.  Presumably  it 
can  be  attributed  to  other  processes  neglected  by  the  model. 


Depth  (m)  Heiaht  (m) 


1  1 1 


Figure  5. 18  The  Model-Predicted  Temperature  Field  at  the 

Rim  for  Experiment  11 
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.  19  Ihe  Model-Predicted  Temperature  Field  over  the 
Slope  with  the  Inclusion  of  the  Neutral  Surface 
Layer  for  Experiment  11 


Figure  5 


TEMPERATURE  CO 


1 13 


Figure  5*20  Progression  of  the  1  meter  Temperature  with 

Time  for  Experiment  1 1 
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CHAPTER  6 


SUMMARY  AND  CONCLUSIONS 


6.1  Summary 

Thirteen  experiments  were  carried  out  in  the  North 
Saskatchewan  River  Valley  in  Edmonton  during  1977  and  1978. 
The  aim  of  these  experiments  was  to  determine  the  major 
characteristics  of  the  valley  microclimate,  and  to  use  this 
information  to  create  a  general  simulation  model  ultimately 
capable  of  predicting  the  pollutant  distribution  in  the 
valley.  The  evolution  of  the  valley  temperature  field  is  an 
important  component  in  such  a  mathematical  model.  This 
thesis  has  concentrated  on  providing  the  temperature 
component  by  developing  a  simple  radiative-conductive  model. 

The  present  radiative-conductive  model  met  with  limited 
success.  However,  the  overall  or  gross  characteristics  of 
the  valley  temperature  field  were  simulated  by  the  model. 
Well-mixed  afternoon  conditions  and  poorly-mixed  evening 
conditions  were  replicated.  Vertical  temperature  gradients 
under  both  stable  and  unstable  conditions  were  reproduced. 


, 
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Stronger  inversions  and  lower  temperatures  over  the  north- 
and  east-facing  slopes  than  over  the  opposing  slopes  were 
accounted  for.  The  timing  of  maximum  temperatures  and  the 
onset  of  cooling  was  good.  The  performance  of  the  model  in 
predicting  a  significant  difference  between  an  urban  and  a 
rural  station  was  disappointing.  However,  the  model  did 
predict  slope-rim  temperature  differences  which  were  close 
tc  those  obvserved.  In  this  respect,  the  inclusion  of  a 
neutral  layer  above  the  slope  was  partially  successful  in 
accounting  for  the  effect  of  the  slope  wind.  The  viability 
of  a  double-structure  inversion  above  the  valley  in  the 
evening  was  demonstrated.  Generally,  the  model  reproduced 
all  the  main  micr ome teorological  features,  but  did  so  with 
limited  accuracy. 


6 . 2  Conclusions 

The  most  important  conclusion  concerns  the  application 
of  flat-plane  profile  theory  to  an  inclined  surface.  It  has 
been  demonstrated  that,  for  gentle  slopes,  the  use  of  the 
flat-plane  theory  is  quite  adequate  in  the  well-mixed 
unstable  regime.  In  this  situation,  the  slope  wind  was 
weak,  and  the  error  in  neglecting  the  advective  component  is 
small.  However,  during  the  evening  when  the  boundary  layer 
becomes  stable,  and  cooler  draining  air  initiates  the  slope 
wind,  the  use  of  flat-plane  theory  becomes  questionable. 

Once  the  slope  wind  becomes  well  established,  the 
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logarithmic  wind  profile  of  the  flat-plane  theory  is  not 
observed  in  the  lowest  levels.  In  an  attempt  to  correct 
this  discrepancy,  a  second  wind  profile  was  introduced  in 
the  domain  of  the  slope  wind.  This  profile  was  deemed  to 
have  properties  similar  to  the  flat-plane  profile  with  the 
only  difference  being  one  of  scale.  Due  to  the  mechanical 
nature  of  the  motion,  the  profile  was  assigned  neutral 
stability  characteristics.  The  results  have  indicated  that 
this  procedure  was  reasonably  successful. 

In  addition,  the  importance  of  advection  and  drainage 
flow  has  been  inferred  from  the  results.  The  predicted  late 
evening  slope  temperatures  were  consistently  higher  than  the 
observed.  It  is  felt  that  this  difference  was  mainly  due  to 
the  absence  of  an  advective  component  in  the  model,  and  that 
the  inclusion  of  an  advective  term  would  improve  the 
results. 

The  model  radiative  components  were  generally  adequate. 
On  fair-weather  days,  the  atmospheric  long-wave  flux 
densities  calculated  from  morning  and  afternoon  radiosonde 
observations  differed  by  less  than  2  per  cent.  This  is  less 
than  the  accuracy  of  the  computation.  Hence  the  assumption 
of  a  constant  atmospheric  flux  density  is  good.  The  effects 
of  excluding  specific  humidity  as  a  predicted  variable  in 
the  model,  and  using  a  rough  parameterization  for  the  latent 
heat  flux  density  are  not  completely  known.  Although  the 
flux  densities  were  not  unreasonable,  this  remains  an  area 
for  future  improvement.  On  ideal  cloudless  days  the  solar 
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radiative  flux  density  term  was  adequate.  The  time  and  rate 
of  extinction  of  the  solar  flux  density  was  important.  In 
this  respect,  the  model  was  sensitive  to  the  inclination  and 
orientation  of  the  valley  slope.  Otherwise  the  model  was 
not  sensitive  to  the  slope  parameter.  It  is  conjectured 
that  this  conclusion  would  be  invalid  if  an  advective  term 
was  included.  The  soil  flux  density  exhibited  little 
sensitivity  to  the  value  chosen  for  the  soil  diffusivity. 
This  was  heartening  since  the  soil  parameter  was  assumed 
constant  with  depth  and  independent  of  soil  moisture. 

The  model  proved  sensitive  to  wind.  Whenever  an 
observed  wind  was  input  into  the  model  (every  30  minutes) , 
it  was  reflected  in  the  output-  It  is  suspected  that  this 
sensitivity  was  mostly  due  to  model  shock,  and  that  the 
incorporation  of  wind  as  a  predicted  variable  would 
alleviate  the  problem. 


6 . 3  Scope  for  Future  Work 

In  any  work  involving  boundary-layer  meteorology  the 
inter-relationshi ps  between  the  temperature  and  wind  are 
complex.  Thus  the  determination  of  stability  functions, 
such  as  the  Bichardson  number,  should  not  be  based  on  the 
explicit  separation  of  these  two  variables.  In  this 
respect,  a  priority  for  future  work  is  the  incorporation 
into  the  model  of  wind  as  a  prediction  element.  Moreover, 
the  model  could  be  expanded  into  two  dimensions  to  include 
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the  effects  of  advection,  and  the  drainage  and  pooling  of 
cooler  air.  Another  possibility  for  future  work  is  to 
upgrade  the  radiative  components  with  the  inclusion  of  both 
cloud  and  specific  humidity.  Inclusion  of  cloud  would  be 
important  not  only  in  adjusting  the  short-  and  long-wave 
radiation  received  at  the  surface  but  also,  in  the  case  of 
cumulus  development,  as  an  indicator  for  refining  the  latent 
heat  flux  density  corrections.  Also  important  are  the 
extensions  to  a  polluted  atmosphere  in  which  a 
radiative-flux  divergence  term  would  be  necessary;  and  to 
winter  inversion  conditions.  The  incorporation  of  wind  as  a 
predicted  variable  is  likely  to  be  the  most  important  of  all 
the  suggestions.  Its  inclusion  has  the  potential  to  improve 
dramatically  the  accuracy  in  the  performance  of  the  model. 
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APPENDIX  A 

GRAPHICAL  INTEGRATION 
OF  THE  ATMOSPHERIC  LONG-WAVE  FLUX  DENSITY 


In  section  (2. 4)  the  atmospheric  long-wave  flax  density 
due  to  water  vapour  that  reached  the  earth’s  surface  was 
expressed  by 


LW 


«4  ©  If 


(A.1) 


where  E  is  the  emissivity  due  to  water  vapour  and  w  is  the 
corrected  optical  depth.  Strictly  speaking,  this  equation 
is  only  valid  in  a  plane  stratified  atmosphere.  It  is 
derived  (see,  for  example,  Fleagle  and  Businger,  1963)  by 
summing  the  monochromatic  flux  density  contributions  from 
all  the  differential  volume  elements  at  an  arbritrary  level 
over  all  levels  and  emitting  wave  lengths.  Using  the 
hydrostatic  equation  in  conjunction  with  the  definitions  of 
specific  humidity  and  differential  optical  thickness,  this 
equation  may  be  recast  as 


wv 


a 
S  J 


po  ,  ,3E 


o 


^4© 


(S.2) 


The  rate  of  change  of  flux  emissivity  with  optical  thickness 
is  experimentally  known,  and  the  vertical  profiles  of 
specific  humidity  and  temperature  as  functions  of  pressure 
can  be  obtained  from  an  atmospheric  sounding. 


Thus  equation 
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(A. 2)  may  be  computed  numerically. 

Alternatively,  the  quantity  Q  =  —  qT4  (|^-)  may  be  plotted 
as  a  function  of  the  pressure,  and  the  area  under  the  curve, 
which  represents  the  flux  density,  determined  graphically. 
The  Q— P  plot  is  not  unique.  By  a  suitable  change  of 
variables,  Elsasser  (1942)  constructed  a  nomogram,  known  as 
the  Elsasser  diagram,  in  which  the  abscissa  was  a  function 
of  temperature,  and  the  ordinate  a  function  of  the  corrected 
optical  depth.  This  chart  was  particularly  convenient:  the 
isotherms  were  straight  vertical  lines,  as  were  the 
isopleths  of  zero  and  infinite  corrected  optical  thickness, 
further  the  isopleth  w=0  coincided  with  the  upper  edge  of 
the  chart,  and  the  isopleth  w=~  with  the  horizontal  axis. 
Standard  curves  of  constant  optical  thickness  are  pre-arawn 
on  this  chart.  A  schematic  version  of  the  Elsasser  diagram 
is  given  in  Figure  A.1.  Properties  of  this  diagram  have 
been  detailed  by  Elsasser.  The  only  property  pertinent  to 
present  purposes  is  that  the  cross- hatched  area  in 
Figure  A.  1  represents  the  downward  atmospheric  long-wave 
flux  density  due  to  water  vapour.  Recall  allowance  has 
previously  been  made  (Section  2. 4)  for  the  long-wave  flux 
density  due  to  carbon  dioxide. 

In  determining  the  corrected  optical  depth,  a  square 
root  pressure  correction  was  used.  Newer  empirical  evidence 
indicates  that  a  linear  correction  may  be  more  appropriate. 
However,  Houghton  (1954)  states  that,  in  the  computation  of 
the  net  flux  density,  the  methods  differ  by  only  about  3  per 
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Corrected  Optical  Depth  (g/cm2) 
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Temperature  (°C) 


Figure  A.1  The  Elsasser  Diagram 
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cent.  Thus  Elsasser*s  original  correction  was  used. 

For  model  purposes,  the  atmospheric  radiation  was 
assumed  constant  over  the  forecast  period,  and  was  computed 
using  the  radiosonde  data  taken  from  Stony  Plain.  The 
specific  humidity  was  graphically  calculated  from  the 
temperature  and  dew  point  data  using  a  tephigram. 

Tables  A.  1,  A. 2,  and  A. 3  list  the  radiosonde  data  and  the 
corresponding  calculated  specific  humidity  and  optical  depth 
for  Experiments  8,  9,  and  11.  Using  these  tables,  the  w-T 
curves  were  plotted  on  an  Elsasser  diagram,  and  the 
long-wave  flux  density  determined  with  a  planimeter. 
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Table  A- 1  The  Atmospheric  Sounding  for  Experiment  8 


Pressure 

(mb) 

Temp. 

(°C) 

Dew 

Point 

(°C) 

Specific 

Humidity 

(g/Kg) 

Corrected 

Optical  Depth 
(g/cm2) 

9  22 

24.8 

1  0.8 

10.  4 

0.00 

870 

18.6 

6.  6 

8.0 

0.  46 

712 

4.0 

-  0.3 

5.5 

1.43 

655 

-  0.5 

-  9.5 

3.2 

1.64 

629 

-  3.1 

-23. 1 

1.3 

1.69 

613 

-  2.9 

-22.  9 

1.4 

1.70 

493 

-13.5 

-32.  5 

0.  70 

1.80 

479 

-13.9 

-27.  9 

0.  98 

1. 81 

450 

-17.9 

-35.9 

0.  53 

1.82 

400 

-25.5 

-34.  5 

0.59 

1.84 

313 

-39.3 

-47.  3 

0. 13 

1.86 

Table  A. 2 

The  Atmospheric 

Sounding  for 

■  Experiment  9 

Pressure 

Temp. 

Dew 

Specific 

Corrected 

(mb) 

(°C) 

P  oint 

Humidity 

Optical  Depth 

(°C) 

(g/Kg) 

(g/cm2) 

922 

24.6 

8.  4 

9.4 

0.00 

808 

14.0 

3.0 

6.8 

0.  88 

697 

3.4 

-  1.6 

5.  2 

1.46 

680 

2.0 

-  9.0 

3.3 

1. 53 

577 

-  7.1 

-1  8.  1 

1.9 

1.74 

462 

-16.1 

-34.  1 

0.  56 

1.85 

305 

-40.  1 

-54.  1 

0.06 

1. 88 

' 

X 

• 

127 


lable  A. 3  The  Atmospheric  Sounding  for  Experiment  11 


Pressure 

(mb) 

Temp. 

<°C) 

Dew 

Point 

<°C) 

Specific 

Humidity 

(g/Kg) 

Corrected 

Optical  Depth 
(g/cm2) 

930 

1  9.6 

6.6 

7.7 

0.  00 

850 

11.8 

2.8 

6.  1 

0.  53 

828 

9.6 

1.  6 

5.8 

0.65 

778 

5.8 

-5.8 

3.  8 

0.87 

700 

-  0.1 

-  9.  1 

3.  1 

1.11 

637 

-  5.7 

-13.  7 

2.4 

1.25 

589 

-  8.9 

-27.  8 

0.84 

1.31 

574 

-  9.9 

-20.  8 

1.5 

1.32 

500 

-17.1 

-32.  1 

0.62 

1.38 

400 

-28.7 

-44.6 

0.  15 

1.41 

376 

-32.  3 

-36.  1 

0.50 

1.42 

321 

-40.1 

-44.6 

0.  15 

1.43 

■ 
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APPENDIX  B 

THE  SOLAR  FLUX  DENSITY 

In  section  (2.4)  the  expression  used  to  evaluate  the 

solar  flux  density  was  given  as 

F  =  Sn(l-jj){  sin(6)cos(£)  +  cos  (6)  sin(£)  cos  (HAN-wt-HATO)  } 

SW  u 

(B.1) 

where  SQ  =  solar  constant 

v  =  albedo-turbidity-diffuse  solar  radiation  factor 
<5  =  declination  of  sun 

£  =  angle  between  celestial  pole  and  slope  normal 

HATO  =  hour  angle  of  sun  at  t  =  0 
HAN  =  hour  angle  of  slope  normal 
w  =  earth* s  angular  rotation 
t  =  time 

This  eguation  results  from  the  cosine  law  for  radiation, 
na  mely 

F  =  Fq  cos (y)  (B .  2) 

where  F0  is  the  radiative  flux  density  through  a  surface 
when  the  flux  beam  is  normal  to  that  surface,  and  Y  is  the 
angle  of  incidence  between  the  beam  and  the  surface  normal. 
It  is  clear  that  F0  can  be  identified  with  S0*(1-u).  Thus 
it  remains  to  show  that  the  angle  of  incidence  is  given  by 

y  =  cos  ^  {  sin(6)cos(£)  4-  cos (6) sin (£) cos (HAN-wt-HATO)  } 

From  the  spherical  geometry  illustrated  in  Figure  B.1,  the 
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ZE 


CP  =  celestial  pole 


ZE  =  zenith 


H  =  hour  angle  of  slope 
normal  relative  to 
hour  angle  of  sun 


a  =  slope  inclination 
3  =  slope  orientation 


HAN  =  hour  angle  of  slope 
normal 

HAS  =  hour  angle  of  sun 

N  =  north 

S  =  sun 

SN  =  slope  normal 


y  =  angle  of  incidence 

between  solar  flux 
beam  and  slope 

6  =  declination  of  sun 

A  =  latitude 

£  =  angle  between  celestial 

pole  and  slope  normal 


ZD  =  zenith  distance 


Z  =  south  point 


Figure  B.  1  The  Celestial  Dome 


following  relationships  hold 
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cos(y)  =  sin(6)cos(0  +  cos  (6)  sin(  £)  cos  (H) 


(B.4) 


H 


HAN  -  HAS 


(B.5) 


If  the  model  is  started  at  time  t  =  0,  then  the  hour  angle 
of  the  sun  will  be 


HAS 


(B .  6 ) 


wt  +  HATO 


where  HATO  is  the  hour  angle  of  the  sun  at  t  =  0  (initial 
hour  angle).  Wnen  equations  (B.5)  and  (B.6)  are  combined 
with  equation  (B. 4) ,  the  expression  for  the  angle  of 
incidence  given  by  eguation  (B.3)  results. 

In  the  model,  the  zero  hour  was  chosen  as  1235  MST. 

This  time  is  very  close  to  local  noon,  and  hence  the  initial 
hour  angle  is  very  close  to  zero.  Values  for  both  the  sun's 
declination  and  initial  hour  angle  were  extracted  from  the 
1977  Nautical  Almanac10  and  corrected  appropriately. 

Sunrise  or  sunset,  neglecting  both  atmospheric 
refraction  and  the  semi-diameter  of  the  sun,  will  coincide 
with  a  solar  flux  density  of  zero*  From  eguation  (B.4), 
this  must  occur  at 


+  HAN 


(B.7) 


HAS 


This  expression,  in  conjunction  with  eguation  of  time  and 
longitude  corrections,  was  used  to  determine  sunrise  and 
sunset  over  the  slope.  Values  for  the  equation  of  time 


10published  by  the  U.S.  Government  Printing  Office, 
Washington,  D.C. 
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corrections  were  also  taken  from  the  almanac.  Between 
sunset  and  sunrise,  the  calculated  solar  flux  density  was 
reset  to  zero. 
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APPENDIX  C 


GAUSSIAN  ELIMINATION  OF  A  TBIDIAGONAL  MATEIX 


A.Y  +  B.Y.  +  C.Y.  =  D.  ,  j=  1 ,  .  .  .  ,N 

3  3+1  J  3  3  J-l  3 


This  appendix  describes  the  solution  by  Gaussian 
elimination  of  the  set  of  equations 

(C.  1) 

where  and  Y^  are  known. 

Since  YN  is  known,  equation  (C.1)  can  be  solved  for 

j  =  N-1  in  terms  of  one  unknown  Y^-2  •  ie 

Y  =  F  +  G  Y 
N-1  N-1  N-1  N-2  ’ 


where  F 


N-1 


PN-1  \-l  fn 

B”  *  +  V’  G'  ’ 


-  C 


N-1 


N-1 


1  N 


N-1 


B”  ’  +  V:  a 


N-1 


•1  N 


F  =  Y 
N  N  ’ 


and 


gn  -  0 


Using  this  expression  for  YN_^  ,  equation  (C. 1)  for  j=N-2 
becomes 

Y, 


'N-2  FN-2  +  GN-2  YN-3 


D 


where  F 


N-2 


F 


2  N-1 


-  C 


N-2 


N-2 


B~  -  +  V-  G 


N-2 


2  N-1 


N-2 


B"  -  +  V-  G 


N-2 


2  N-1 


This  process  of  determining  Y.j  in  terms  of  one  unknown  Y^_1 
is  continued  yielding 


Y.  =  F.  +  G  Y  , 

3  J  J  J"1 


where  f 


D.  -  A.  F. 

3  3 J+L 

B.  +  A.  G 
3  3  3+1 


-  C 


Gj  Bj  +  k..  G. 


3  3  3  +  1 

Ultimately  j  =  2  is  reached.  Since  Y-  1=1  is  known,  Y2  is 
determined.  Back  substitution  of  Y2  then  yields  Y3.  This 


. 
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process  of  back  substitution  is  continued  until  all  values 

of  y.  are  determined. 

3 


134 


APPENDIX  D 

THE  MODEL  GRID 

The  values  for  the  grid  points  employed  in  all 
computer  simulations  are  given  in  Tables  D. 1  and  D. 2. 


Table  D. 1  The  Grid  Point  Values  within  the  Boundary  Layer 


Grid 

Point 

Height 

(m) 

Grid 

Point 

Height 

(m) 

1 

0.  0 

18 

5.  000 

*  10  +  2 

2 

1.000 

★  IQ*2 

i  9 

6.  000 

*10+2 

3 

2.  1  54 

*10-2 

20 

7.  000 

*1  0  +  2 

4 

4.  642 

*‘l  Q-2 

21 

8.  000 

*1  0  +  2 

5 

1.000 

*10-1 

22 

9.  000 

*10  +  2 

6 

2.  154 

*oo-i 

23 

1. 000 

*  10  +  3 

7 

4.642 

*10-1 

24 

1.  100 

*10+3 

8 

1.  000 

25 

1.  200 

*10+3 

9 

2.  154 

26 

1. 300 

*10  +  3 

10 

4.  642 

27 

1.  400 

*  10  +  3 

11 

1 .000 

*1  0*1 

28 

1. 500 

*1  0+3 

12 

2.  1  54 

*'!  0+i 

29 

1.  600 

*1  0+3 

13 

4.  642 

*0  0+  1 

30 

1.  700 

*10  +  3 

14 

1.000 

*10  +  2 

31 

1.  800 

*10  +  3 

15 

2.000 

*1  0  +  2 

32 

1.900 

*1  0+3 

16 

3.  000 

*1  0  +  2 

33 

2.  000 

*10+3 

17 

4.  000 

*10  +  2 
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Table  D.2 

The  Grid  Point  Values 

within 

the  Soil 

Layer 

Grid 

Point 

Depth 

(m) 

Grid 

Point 

Depth 

(m) 

1 

0.  0 

5 

-4.472 

*10-2 

2 

-5.  0  0  0  *1  0~3 

6 

-9. 283 

*10-2 

3 

-1.  038  *10-2 

7 

-1.927 

*10-i 

4 

-2. 154  *10-2 

8 

-4.000 

*10—1 

- 
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APPENDIX  E 

MODEL  PARAMETERS 


Table  E. 1  details  the  model  parameters  common  to  all 
computer  simulations. 


Table  E. 1  Common  Model  Parameters 


Parameter 

Programme 

Name 

Value 

Latitude 

ALAT 

53°  33 ' N 

Longitude 

ALON 

113°  30 ' W 

Slope  Elevation 

A 

16.25° 

Slope  Orientation 

B 

103°  E 

Reference  Height 

ZNO 

0.  01  m 

Height  of  Sub-Layer 

ZT1 

100.  m 

Height  of  Boudary  Layer 

ZT2 

2000.  m 

Number  of  Points  between 

and  including 

a)  surface  and  ZNO 

NPGNO 

2 

b)  ZNO  and  ZT 1 

NPN0T1 

13 

Grid  Spacing  above  Sub-Layer 

DZ2 

100  m 

Depth  to  Soil  Bottom 

ZSB 

-0.4  m 

Soil  Reference  Depth 

ZNOS 

—0.005  m 

Number  of  Points  between 

Surface  and  Soil  Bottom 

NZSPTS 

8 

Start-up  Time 

TO 

12  hr  35  min  MDT 

Termination  Time 

TF 

24  hr  35  min  MDT 

Time  Step 

D1 

5  min 

Solar  Constant 

SO 

1353  W  m-2 

Carbon-dioxide  Correction 

CF 

.  82 

Albedo-Turbidity-Dif f use 

Solar  Radaition  Factor 

ALB 

.39 

Stephan-Boltzmann  Constant 

SIGMA 

5.67*10-8  Wm— 2 °K~4 

Earth's  Angular  Rotation 

W 

7.292*10-5  sec* 1 

. 


137 


Table  E. 1  Continued 


Parameter 

Programme 

Name 

Value 

Eddy  Viscosity 

KMOL 

2.2*10-5  m-ssec-1 

Soil  Diffusivity 

KHS 

1.5*10~7  m-2sec_1 

Density  *  Heat  Capacity 

LAM 

1.2*10+3  Jm-30K-i 

Soil  Conductivity 

LAMS 

2.5*10-1  Wm-i°K-i 

Error  Limit  for  Surface 

DER2 

.00001  °C 

Temperature  Iteration 

Maximum  Number  of  Iterations 

NITB 

100 

for  Surface  Temperature 

Error  Limit  for  Temperature 

DEE  4 

.00001  °c 

Profiles  Iteration 

Maximum  Number  of  Iterations 

NITC 

500 

for  Profile  Determination 

Weights  for  Implicit  Method 

a)  for  Air 

A  L, 

BE 

1.5,  0.5 

b)  for  Soil 

ALS, 

BES 

Displacement  Time  for 

TM 

.  5  hrs 

Friction  Velocity  Data 

Time  Interval  between 

TL 

.  5  hrs 

Friction  Velocity  Data 

Von  Karman’s  constant 

VKC 

0.  4 

' 
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APPENDIX  F 

THE  COMPUTER  PROGRAMME 

The  following  pages  of  this  appendix  list  the  FORTRAN 
code  for  the  computer  model. 


. 


Co  Qn  Co  CO  Co  co 
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IMPLICIT  RE AL*8 ( A-H,0-Z) ,  INTEGER (I- N) 

LO  GIC A  L* 1  F  (1) /' **/ 

DOUBLE  PRECISION  KN  (10  )  ,  KNN  ( 50  ,  1  0)  ,  KH  ( 1  0)  ,  K  HH  (  50 , 1 0)  , 
&  LAM,LAMS,KMOL,KHS 

COMMON  /A  1/1, J 

&  /A2/TO,TF,DT,T,TM,TL 

6  /A3/ALAT, ALON 

&  /A4/A(1 0) ,B(10) , IDATE,IDUM1 

&  /A5/SST  (10)  ,  SET  (10)  rZET(10)  ,  HAN  (10)  ,  DEC,  EOT 

&  /A6/NSTNS  ,NZPTS ,  NZPTS1  ,NZ  PTSA,  NPGN0  ,NPN0T  1, 

6  NZSPT  S,  N  ZSPT  1 

&  /A7/ZN0,ZT1 ,ZT2,DZ2,ZN0S,ZSB 

&  /A8/Z  (50)  ,ZS  (10) 

8  /A9/THET (2,50, 10) ,TS (2,10,10) 

8  /A10/H22, H12,HX, HS22, HS 1 2 , HSX 

&  /A  1  1/LAM, LAMS, SIGMA ,CF , DER2 ,NITB, NIT2 

&  /A  1 2/KN, K NN  /A  1 3/KH, KHH 

&  /A14/FLWWV (10) ,FSW,FSH,FSHS,FLH,FLW 

&  /A1 5/DZ  A ( 50)  ,DZB  (50) ,DZC (50) ,C2 (50) ,C3 (50) 

&  /A 16/DZ  AS  (10)  ,DZBS(10)  ,DZCS  (10) ,C2S  (10)  ,C3S(10) 

&  /A  17/SO, W,HATO,ALB, DER4 , NITD, NIT4 

&  /A18/USTAR(10)  ,V  KC , KMO  L 

&  /A  1 9/RI , PHI, BRI , UB 

&  /A20/K, IDUM 

&  /A 2 1/PHILIM, TSLP, JSLP,  IDUM2 

C 
C 

DIR  (X)=X*.  17453292519 94 329 D- 01 
EDI  (X) =X/.  1745329251994329D-01 

D2R  (X)  =  (AINT  (SNGL  (X)  )  +( X-AINT  (SNGL  (X)  )  )  *  1 00.  D0/60.  DO) 

&  *. 1745329251994329D-01 

D3R  (X)  =  (  AINT  (SNGL  (X)  ) 

&  +(  (X-AINT  (SNGL(X))  )  *  1 00, D  0 

-  (X*1 00, DO- AINT (SNGL (X)  *100,)  )  ) /6  0.D0 
♦  (  X*100.  DO -AINT  (SNGL  (X)  *100,)  )/36.D0  ) 

*,  17  4532925  I  9943 29D-01 

H2H  (X) =AINT (SNGL  (X) ) + (X-AINT (SNGL (X) ) ) *  100. D0/60 .DO 
H3H  (X)  =  AINT  (SNGL  (X)  ) 

+  (  (X-AINT  (SNGL  (X)  ))  *100.  DO 

-  (X* 1 0  0. DO- AINT (SNGL (Xj *100.  )  )  ) /60. DO 
+(  X*100. DO-AINT(SNGL (X) *100, )  ) /36. DO 

C 
C 

CALL  SUN 
C 
C 

DO  100  1  =  1 , NSTN  S 
WRITE  (7, F)  SET (I )  , SST (I) 

100  CONTINUE 
C 
c 

CALL  GRD 
C 


- 


on  noon  no  no  on  on 
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C 


READ  (4, F)  H  ATO,  TO, TF, DT , 


£ 

6 

£ 

£ 

£ 

£ 


S0,W  ,  SIGMA,  CF,  AL  B#  FLWA, 
KMOL,KHS,  LAM,  LAMS, 

DER2 ,DER4 , 

NITB, NITD, 

AL,  BE,  ALS,BES, 
TM,TL,VKC,PHILIM,  TSLP 


T  =0. DO 


1 10 

120 

130 


140 


HATO  =  D2R  (HATO) +D1R  (EOT*15,DO) 

TO  =H2H (TO) 

TF  =H2H  (TF) 

DT  = H2H (DT) 

NSTPS=IFIX ( SNGL  (  (TF-TO)/DT+.  5)  )  +1 


DO  110  I=1,NSTNS 

FLWWV  (I)  =FLWA*DCOS  (A  (I)/2-  DO)  *DCOS (A  (I)/2.D0) 
CONTINUE 


DO  120  J=1 , N  ZPTS 
WRIT  E  ( 7,  F)  Z(J) 
CONTINUE 

DO  130  J= 1,NZSPTS 
WRI TE ( 7, F)  ZS(J) 
CONTINUE 


H22  =  Z (3) *Z(3) 

HI  2  =  Z  (2)  *Z  (2) 

EX  =  (Z  (3)  -Z  (2)  )  *Z  (3)  *Z  (2) 
H522=  Z S (3) *  ZS (3 ) 

HS  12=  ZS  (2)  *ZS(2) 

KSX  = (ZS (3)-ZS (2)) *ZS  (3) *ZS(2) 


CALL  FVEL 


DO  140  J=2 , N  ZPTS 1 
DZA  (J)  =Z(J)  -Z(J-1) 

DZB  (J)  =Z(J+1)-Z(J) 

DZC(J)  =  {DZA  (J) +DZB  (J)  )/2.D0 
C2(J)  =AL*DT*3600-  DO/ 2-  D0/DZC  ( J) 

C3  (J)  =BE*DT*36 00 .D0/2-  D0/DZC  (J) 

CONTINUE 


DO  150  J=2 , NZS PT 1 
DZAS(J)  =ZS(J)  -ZS(J-I) 
DZBS(J)  =ZS  (J+1) -ZS  (J) 


1 10 

120 

130 


140 


. 


no  non  n  non  oonnononn 
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DZCS(J)  =  (DZAS (J)  +  DZBS (J)  ) /2. DO 
C2S  (J)  =ALS*DT*3600. D0*KHS/2 . DO/DZCS (J) 

C3S  (J)  =BES*DT*3600.  D0*KHS/2. DO/DZCS  (J) 

150  CONTINUE 


WRITE (7,F)  NSTNS,N  ZPTS,NSTPS, NZSPTS 


CALL  SSTS 


DO  200  I=1,NSTNS 
CALL  INIT 

FS W=SO*  (  DSIN  (DEC)  *DCCS  (ZET  (I)  ) 

&  +  DCOS  (DEC)  *  DSIN  (ZET  (I)  )  *DCOS  (HAN  (I)  -HA  TO)  ) 

5  * (1- DO- ALB) 

IF( ( (T  +  TO)  . LT.  S ST  (I)  )  . AND. 

6  (  (T  +  TO)  -  GT.  SRT  (I)  )  ;  GOTO  210 

IF  (  (  (T+TO)  .  LT.  (S ST  (I)  +24.  DO)  )  .AND. 

S  (  (T  +  TO)  .  GT.  (SRT(I)  +24.D0)  )  )  GOTO  210 
FSW=0. DO 
210  CONTINUE 

CALL  DFVTY 

CALL  TSFC 

TS  (2,1,I)=THET  (2,  1,1) 

TS  (1,1,1) =THET  (2,1,1) 

THET  (1,1,1) =THET (2,1,1) 

K0=  1 

WRI TE  (7  ,  F)  NIT 2,  I,  KO 

WRITE  (6)  (THET  (2  ,  J,I)  , J=1,NZPTS)  , 

S  (TS  (2  ,  J , I)  ,J=1, NZSPTS) 

FNET=FSW  +  FLWWV  (I)-FLW 

WRI TE ( 8)  FSW,FSH,FSHS,FLH,FLWWV (I)  , FL W, FNET 
200  CONTINUE 


DO  900  K=2 ,NSTPS 

IF (  DMOD (T  +  TM+TO+1 .D-06,TL)  .  LT.  1. D-06  ) CA LL  FVEL 
T=T  +DT 


DC  600  I=1,NSTNS 

CALL  DFVTY 
CALL  DFVTY2 


no  no 
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C 

c 


c 


220 


230 

C 

600 


CALL  PBOF 

WHITE (7, F)  NIT2,NIT4 ,I,K 

WHITE  (6)  (T  RET  (2,  J,I)  , J=1 ,NZPTS)  , 

&  (TS  (2,  J,I)  ,J=1,NZSPTS) 

FNET=FSW+FLWWV  (I)  -FLW 

WRITE  (8)  FSW,FSHrFSHSrFLH,FLWWV (I)  ,FLW,FNET 

DO  220  J=  1  r NZPTS 

THET  ( 1  ,  J  ,  I)  =THET  (2,J,I) 

CONTINUE 

DO  230  J=1 r  NZ  SPTS 
TS  (1  ,J ,1) =TS  (2, J,I) 

CONTINUE 

CONTINUE 


900  CONTINUE 


STOP 

END 
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SUBROUTINE  SUN 

IMPLICIT  REAL*8  (A-H,0-Z) , INTEGER (I-N) 

LOGICAL* 1  F  (1) /■ * •/ 

COMMON  /A1 /I , J 
8  /A3/ALAT,  ALON 

&  /A4/A  (10)  ,B  (10)  ,  IDATE,IDUM1 

&  /A5/SST  (10)  ,SRT  (10)  ,  Z  ET  ( 1  0)  ,  H  AN  ( 1  0)  ,  DEC,  EOT 

&  /A6/NSTNS ,IDUM (7) 

C 

DIR  (X)=X*.  174532 92519 94 329 D- 01 
RD1  (X) =X/.  1 745  32  92519  94329D-0 1 

D2R  (X)  =  (AINT  (SNGL  (X)  )  +( X-AINT  (SNGL  (X)  )  )  *  100,  D0/6  0.  DO) 

5  *.  1745329251994329D-01 
D3R  (X)  =  (  AINT  (SNGL  (X)  ) 

6  +{  (X-AINT  (SNGL  (X)  )  )  *100.  DO 

5  - (X*1 00- D0-AINT(SNGL (X) *100.) )  )/60.D0 

6  +(  X*100. DO-AINT  (SNGL  (X) *100.)  J/36.D0  ) 

&  *.1745329251994329D-01 

H2H (X) =AINT (SNGL (X) ) + (X-AINT (SNGL (X) ) ) *100. D0/60. DO 
H3H  (X)  =  AINT  (SNGL  (X)  ) 

&  +(  (X-AINT  (SNGL  (X)  ))  *100.  DO 

5  - (X*1 00. DO-AINT(SNGL (X) *100. ) )  J/60.D0 

6  +(  X*100.D0-AINT(SNGL (X) *100. )  ) /36. DO 

C 

c 

READ  (4 , F)  NSTNS 

READ  (4 , F)  ALAT, A  LON , DEC, EOT, I  DATE, 

S  (A  (I)  #B(I)  rI=1r NSTNS) 

C 

c 

N  =  0 

ALNCOR=DMDD  (H2H  (ALON)  r  1  5. DO)  /15.D0 
ALAT=D2R(  ALAT) 

ALON=D2R(  ALON) 

DEC  =  D2R  (DEC) 

EOT  =H3  H  (EOT) 

C 

c 

DO  290  1=1,  NSTNS 
A  (I)  =D1R  (A  (I)  ) 

B  (I)  =D1  R  (B  (I)  ) 

IF  (I.GT.  1)  GOTO  130 

IF  (ALAT.  EQ.  DIR  (90.D0)  .  AND.  DEC.  GT.  0.  DO)  GOTO  110 
IF(ALAT.EQ.  DIR  (90.  DO)  )GCTO  100 
ARG  =  -DT  AN  (ALAT)  *  DTAN  (  DEC) 

IF  ( AR  G.  GE-  1 .  D  0)  GOTO  100 
I  F(  ARG.LT.- 1- DO)  GOTO  110 
T  =  RD 1  (ACOS (ARG)  )  /1 5. DO 
SST  (I) =  T+12.DO+EOT+  ALNCOR 
SRT  (I) =-T  + 1 2. D 0  +  EOT+ ALNCOR 
GOTO  120 

100  WRITE  (5,300) 

STOP 

110  WR ITE  (5,310) 


- 
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SST  (I)  =  50.  DO 
SET  (I)=-50.  DO 
120  CONTINUE 

ZET(I) =D1E(90. DO)-ALAT 
HAN  (I)  =0.  DO 
GOTO  290 
130  CONTINUE 
C 
C 

HAN  (I)=0.  DO 

IF(  ( B  (I)  .  EQ.O.  DO.  AND.  A  (I)  .LE.  (D1  E  (  90 .  DO)  - AI  AT)  ) 

8  .OF.  (B(I)  .EQ.D1E(180.DO))  )  GOTO  160 

IF  (B  (I)  .EQ.O.  DO)  GOTO  150 
IF(  DCOS (ALAT) *DCOS  (A  (I) )  .EQ. 

6  DSIN(ALAT)  *DCOS  (B  (I)  )  *DSIN  (A  (I)  )  )  GOTO  140 

HAN  (I)  =DATAN  (  DSIN(A(I)  )  *DSIN{B  (I)  ) 

8  /(  DCOS  (ALAT) *DCOS  (A (I) ) 

8  -DSIN  (ALAT) *DCOS  (B (I) )  *DSIN  (A  (I)  )  )  ) 

IF  (  (B (I) .GT.O.DO. AND. HAN (I) .GT.O.DO) 

8  -OF.  (B (I) . LT. 0 .DO . AND . HAN (I) . LT. 0. DO)  )  GOTO  160 
HAN  (I)  =D1R  (  180.  DO)  +HAN(I) 

IF  (B  (I)  .GT.O.DO)  GOTO  160 
HAN  (I)=HAN  (I)  -DIE  (3  60.  DO) 

GOTO  160 

140  HAN  (I) =D1B{90. DO) 

IF  (B  (I)  .GT.O.DO)  GOTO  160 
HAN  (I)  =-HA N  (I) 

GOTO  160 

150  HAN  {I)  =D1E  ( 180.  DO) 

160  CONTINUE 

IF  (HAN  (I)  .  EQ.  O.DO.AND.B  (I)  .LE.  DIE  (90.  DO)  )  GOTO  170 
I  F  (  HAN  (I)  -  EQ.O.  DO)  GCTO  180 
IF  (HAN  (I) . EQ. DIE  (1 80. DO) ) GOTO  ISO 

ZET  (I)  =DARS IN  (DSIN  (A(I)  )  *DSIN  (  B  ( I)  )  /DSIN  (HAN  (I)  )  ) 

GOTO  200 

170  ZET  (I)  =D1  R  (  90  .  DO ) -ALAT- A  (I) 

GOTO  200 

180  ZET  (I) =D1R  (90. DO  )  -  ALAT+ A  (I) 

GOTO  200 

190  ZET  (I) =A(I) -DIR (90. DO) +ALAT 
200  CONTINUE 

IF  (ZET  (I)  .EQ.D1R  (180. DO),  AND.  DEC*.  GE.  0 ,  DO)  GOTO  210 

IF  (ZET  (I) . EQ. DIR  (180. DO) ) GOTO  220 

IF  (ZET  (I)  .EQ.O.  DO,  AND.  DEC.  GT.O.DO)  GOTO  220 

IF(ZET  (I). EQ.O. DO)  GOTO  210 

ARG=0. DO 

IF  (ZET  (I)  .EQ.  DIR  (90.  DO)  .OR.  ZET  (I)  .  EQ.  DIE  (-90.  DO)  ) 

8  GOTO  230 

ARG=-DTAN (DEC) /DTAN  (ZET (I)  ) 

IF(ARG.GE. 1. DO) GOTO  210 
IF  (ARG.LT.-1.  DO)  GOTO  220 
GOTO  230 

210  WRITE  (5# 320 ) 

STOP 


- 

' 


n  n 
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220  WRITE  (5,330) 

SST(I) =50. DO 
SET (I) =-50. DO 
GOTO  260 
230  CONTINUE 

H  =  ACOS  (AEG) 

T1  =RD  1  (HAN  (I)  -H)  /1 5.  DO 
T  2=ED 1 (H  +  HAN (I) )  /I  5. DO 
SET  (I)=T1  +  12-DO  +  EOT  +  ALNCOR 
SST(I)  =T2+ 1 2- DO  +EOT+A  LNCOR 
IF  (SET(I)  .GT.O.  DO)  GOTO  240 
SET  (I) =SRT (I) +24. DO 
SST  (I)  =SST  (I)  +24. DO 
240  CONTINUE 

IF  (SRT(I)  .LT.  SST  (1)  )  GOTO  260 

IF  (SST  (I) .LE.  (SRT(1) +24. DO)  )  GOTO  250 

SST  (I) =SST (I) -24. DO 

SRT(I)  =SRT  (I)  -  24.  DO 

GOTO  260 

250  WRITE  (5r34  0)  SRT  (I)  ,  SST  (I) 

STOP 

260  CONTINUE 

IF  (SST(I)  .LT.  SST  (1)  )  GOTO  270 
SST (I) =SST (1) 

270  IF  (SRT  (I).  GT.  SRT  (1)  )  GOTO  280 
SET  (I)  =SRT  (1) 

280  CONTINUE 
290  CONTINUE 


RETURN 

300  FORMAT  (//, 25H  SUN 
&  17H  EXECUTION 

310  FORMAT  {//,25H  SUN 
320  FORMAT (//,24H  SUN 

5  17H  EXECUTION 

330  FORMAT  (//, 23H  SUN 
340  FORMAT (//,26H  SUN 

6  12H  BELOW  PLA; 


NEVER  ABOVE  HORIZON, 
HALTED,//) 

ALWAYS  ABOVE  HORIZON,//) 
NEVER  ABOVE  SLOPE, 
HALTED,//) 

ALWAYS  ABOVE  SLOPE,//) 
ABOVE  SLOPE  ONLY  WHEN 
E, 2F1 0. 5) 


END 


■ 
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SUBROUTINE  SSTS 

IMPLICIT  REAL*8  (A-HrO-Z) , INTEGER (I-N) 
L0GICAL*1  F  (1)  /'  *'  / 

COMMON  /A 2/TO,  TF  ,  DT, DUM  (3) 

&  /A5/SST  (10)  ,SRT  (10)  ,DUM2 (22) 

&  /A6/NSTNS ,IDUM  (7) 

C 

C 

DO  120  1= 1 , NSTNS 

K2=-9  9 

K3=- 9  9 

K4=-9  9 

K5=-99 

IF  (TO.GT.  SRT(I)  )  GOTO  100 
1 1  =  IFIX  (SNG  L  (  DSIGN  (1  .  DO  ,  SRT(I)-TO)  )) 
K 2=  IFIX  (SNGL  (  (SRT  (I )  -TO )  /DT)  )  +11 
K3=K2+I1 
100  CONTINUE 

IF  (TF.  IT.  SST  (I)  )  GOTO  110 
I2=IFIX  (SNGL  (  DSIGN(1.D0,  SST(I)-TO)  )) 
K4=IFIX  (SNGL ( (SST (I) -TO) /DT)  )  +12 
K5=  K4  +  1 2 
C 

110  NRITE(7#F)  I,K2,  K3,K4,K5 
C 

120  CONTINUE 
C 

c 

RETURN 

END 
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SUBROUTINE  GRD 

IMPLICIT  RE AL*8 (A-H,0-Z) r INTEGER (I-N) 
LOGICAL* 1  F  (1) /*  * '/ 

COMMON  /A6/NSTNS  rN  ZPTS , NZPTS 1 #NZPTSA,NPGNO# 

5  NPNOT 1 , NZSPTS, N  ZSPT1 

6  /A7/ZN0, ZT1,ZT2,DZ2,ZN0S,ZSB 

&  /A8/Z  (50)  ,ZS  (10) 

&  /A21/DUM(2) , JSLP ,IDUM2 

C 

c 

READ  (4,F) ZN0,ZT1 rZT2,NPGN0rNPN0T1,DZ2, 

&  ZNOS,ZSB,NZSPTS 

C 

NZPTSA=NPN0T1+NPGN0-1 

NZPTS  =IFIX (SNGL  ( (ZT2-ZT1) /DZ2)  ) +NZPTSA 
NZPTS 1=NZPTS-1 
NZSPT1=NZSPTS- 1 
JSLP  =NPGNQ+7 
C 

c 

Z (1)=0. DO 

Dl  = (ZT1/ZN0) **  (1, DO/DFLOAT (NPN0T1-1)  ) 

ZLI  =DLOG ( D 1 ) 

DO  100  J=2 , NZ PTS A 

Z (J) =ZN0*DEXP (DFLOAT ( J-NPGNO)  *ZLI) 

100  CONTINUE 
C 

NZPTSB=NZPT  S A+ 1 
DO  110  J=NZPTSB, NZPTS 
Z  (J)  =Z  (J-1)  +  DZ2 
110  CONTINUE 
C 

c 

ZS (1) =0. DO 

D2  =  (ZSB/ZNOS)  **  (1,  DO/DFLOAT  (NZSPTS-2)  ) 
ZLSI  =DLOG(D2) 

DO  120  J=2 , NZS PTS 

ZS  (J) =ZN0S*  DE  XP (DFLOAT (J-2)  *ZLSI) 

120  CONTINUE 
C 

RETURN 

END 
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C 

C 


500 


c 

c 


c 


100 

110 


120 

130 

140 


150 


160 

170 


180 

190 

200 


SUBROUTINE  INIT 

IMPLICIT  REAL*8 (A-H,0-Z) , INTEGER (I-N) 

LOGICAL* 1  F ( 1 ) /*  *' / 

DOUBLE  PRECISION  ZI  (50)  , TH (50 )  ,ZIS(10)  ,THS  ( 10) 
COMMON  /A1/I,J 

&  /A6/NSTNS,NZPTS,NZPTS1 ,NZPTSA ,NPGN0 , 

&  NPN0T1,NZSPTS,NZSPT1 

&  /A8/Z  (50)  ,ZS  (10) 

&  /A9/THET { 2#  50,10) ,  TS (2,10,10) 


R  EA D  (4  ,  F)  FLAG7 
IF (FLAG7.NE.1. DO) GOTO  500 

READ  (2)  (THET ( 1 , J,I)  ,J=1,NZPTS)  , (TS (1,J,I)  ,J=1,NZSPTS) 
READ  (2)  (THET  (2, J,I)  ,  J=1,NZPTS)  ,  (TS(2,  J,I)  ,  J=  1  ,  NZSPTS) 
RETURN 
CONTINUE 
J  2=N  ZPTS 

READ  (4  ,  F)  N IDP  ,  (ZI  (J)  ,18(1)  ,  J=1,NIDP)  , 

&  NIDPS,  (ZIS  (J)  ,  THS (J) , J=1, NIDPS) 


DO  200  JI=2,NIDP 

GAM=(TH(JI)  -TH(JI-I)  )  /  ( ZI  ( JI- 1) -ZI  ( JI )  ) 

DO  100  J3=1,J2 
J  4=J3 

IF  (ZI  (  JI)  .LE.  Z  (J3)  )  GOTO  110 
CONTINUE 

IF  (JI.  GT.  2)  GOTO  140 
DO  120  J5=J4,J2 
J6=J5 

IF  (ZI  (JI-1)  .EQ.Z  (J5)  )  GOTO  140 

IF  (ZI  (JI-1)  .LT.Z  (J5)  )  GOTO  130 

CONT INUE 

J  6  =  J  6  - 1 

CONTINUE 

DO  150  J7=J 4, J6 

THET  (1  ,  J7 , I )  =  TH  (JI)  -GAM*  (Z  (J7)  -  ZI  (JI)  ) 
CONTINUE 

IF  (JI  . NE  . 2 )  GOTO  170 
J6A= J6+1 

DO  160  J8  =  J6A,NZPTS 

THET(1  ,J8,I)=TH  (JI)-GAM*  (Z  (J8)  -ZI  (JI)  ) 
CONTINUE 

IF  (JI. NE. NIDP)  GOTO  190 
IF  (ZI  (JI;  .  EQ.Z  (1)  )  GOTO  200 
DO  180  J9=1,J4 

THET  (1,J9,I)=TH(JI)  -GAM*  (Z(J9)-ZI  (JI)  ) 
CONTINUE 
J2=J4 
J6=J  4- 1 
CONTINUE 


C 


' 


' 


o  o  n 
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J2=NZSPTS 

DO  400  JI=2,NIDPS 

G  AMS=  (THS(JI)  -THS  (JI-1)  )/(ZIS  (JI-1)  -ZIS  (JI)  ) 

DO  300  J3=1#J2 
J4=J3 

IF  (ZIS (JI) . GE. ZS  (J3) ) GOTO  310 
300  CONTINUE 

310  IF  (JI.GT. 2) GOTO  34  0 

DO  320  J5=J4,J2 
J6=J5 

IF  (ZIS (JI-1) - EQ. ZS  (J5) ) GOTO  340 
IF  (ZIS  (JI-1),LT.  ZS(J5))  GOTO  330 
320  CONTINUE 

330  J6=J6-1 

340  CONTINUE 

DO  350  J7=J4,J6 

TS  (1,  J7,I)  =THS  (JI)  -GAMS*  (ZS(J7)  -  ZIS  (JI)  ) 
350  CONTINUE 

IF  (JI .  NE.  2)  GOTO  37  0 
J6A=J6+ 1 

DO  360  J 8=J 6A , NZ  SPTS 

TS  (1  ,  J8 ,1)  =THS  (JI)  -GAMS*  (ZS(J8)  -  ZIS  (JI)  ) 
360  CONTINUE 

370  IF  (JI  .NE.NI  DPS)  GOTO  390 

IF  (ZIS  (JI)  .EQ.  ZS(1)  )  GOTO  400 
DO  380  J9= 1 r J4 

TS(1,  J9#I)  =  THS  (JI)  —GAMS*  (ZS  (J9)  -ZIS  (JI)  ) 
380  CONTINUE 

390  J  2= J4 

J6=J4- 1 
400  CONTINUE 


DO  520  J  =  1 f  N  ZPIS 
THET (2, J, I)  =THET  ( 1 , Jr I) 
520  CONTINUE 

DO  530  J=1,NZSPTS 
TS(2,  J,I)  =TS  ( 1  /  J  #  I) 

530  CONTINUE 
C 

RETURN 

END 


. 


' 


no  n  n  no 


SUBROUTINE  TSFC 
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IMPLICIT  R  EAL*8  (A-H,0-Z) , INTEGER (I-N) 
LOGIC AL*1  F ( 1 ) /' **  / 

DOUBLE  PRECISION  KH  (1 0)  ,KHH (50, 1 0)  , 


& 


LAM, LAMS ,  KMOL,  KHS 


COMMON  /A 1 /I , J 


5 

6 
& 
8 
& 
& 
s 


/A2/TO,  DUM4  (2)  ,T,DUM5  (2) 

/A8/Z  (50)  ,ZS  (10) 

/A9/THET  (2,50,10) ,TS (2,10,10) 
/A10/H22, H12,HX,HS22,HS12,HSX 
/A  11/LAM, LAMS, SIGMA,CF ,DER2 ,NITB, NIT2 
/A13/KH, KHH 

/A14/FLWWV  (10)  ,FSW,FSH,FSHS,FLH,FLW 


100 

110 

115 

C 

C 

117 


NIT  2=0 
FL  AG2=0 . DO 
CONTINUE 
NIT2=NIT2+1 

IF (NIT2. LT. NITB) GOTO  110 
T2=TO+T 

WRITE(5,180)  NIT2,T2,I,ERR2 

STOP 

CONTINUE 


FSHS=-LAMS*(  H  S2  2*  TS ( 2 , 2,1) - HS 1 2*TS (2,3,1) 

&  -  (HS22-HS12)  *TS(2,  1,1)  ) /HSX 

FSH=LAM*KH (I) * (  H22*THET  (2,2,1) -HI 2*THET (2 , 3 , I) 

&  -  (H22-H12)  *THET(2,  1,1)  )  /KX 

FL  H=FSH 

IF(FSH.LT.O.DO) GOTO  115 

FLH=-.5D0*FSH 

CONTINUE 

FL  W=FSW+FSH+FSHS+FLWWV  (I)  +FLH 
IF (FLW.LT.0.D0) GOTO  120 

IF  (NIT2.GT. 20  .AND,  DABS  (FSH)  . LT. 20. DO  ) GOTO  200 
FL  A  G2= 1 . DO 

GT=THET (2,1,1) - ( FLW/SIGMA/CF) **. 25D0 
GTP=1.  D0+.  2  5D0*  (  (LAMS/ (-ZS  (2)  )  +2.  DO* LAM*  KH  (I)  /Z  (2)  ) 
8  /SIGMA) /(FLW/SIGMA/CF) ** (.75D0) 

IF (FSH.LT.O. DO) GOTO  117 

GTP=1.D0+.2  5D0*  (  (LAMS/(-ZS  (2)  )  +-  5 DO*  LA M*  KH  (I)  /Z  (2)  ) 
8  /SIGMA) / (FLW/SIGMA/CF) ** (.  75  DO) 

CONTINUE 

THETF=THET (2, 1 ,1) 

ERR  2=  GT/GT  P 

THET(2,1,I)  =  THET  (2,1,I)-ERR2 
IF (DABS (ERR2) . LT. DER2) GOTO  160 
GOTO  100 


100 

110 

115 

C 

C 

117 


151 


C 

120  CONTINUE 

GT2=THET  (2, 1,1) 

IF  (FLAG 2.  EQ-  1 .  D  0)  GOTO  150 
C 

THETF=250. DO 
130  CONTINUE 

FSHS  =  LA  MS/ZS (2) *  (THETF-TS (2 , 2 , 1)  ) 
FSH=LAM*K H  ( I)  /Z  (2)  *  (THET  (2,2,1)  -THETF) 
FLH=FSH 

IF  (FSH. LT.  0.  DO)  GOTO  135 
FLH=-.5D0*FSH 
135  CONTINUE 

FLW=FSW+FSH+FSHS+FL  WWV(I)  +FLH 
IF  (FLW-.GT.  0.D0)  GOTO  140 
I F ( TH  ETF.  EQ .  200.  DO ) GOTO  170 
THE  TF  =2  00.  DO 
GOTO  130 
140  CONTINUE 

GT=THETF- (FLW/SIGMA/CF) **.25D0 
C 

150  CONTINUE 

ERR  2=  (GT2-THETF)  /(GT2-GT)  *GT2 
THET  (2, 1 ,1)  =GT2-ERR2 
IF  (DABS  (ERR2)  .LT,DER2) GOTO  160 
GOTO  100 
C 

200  CONTINUE 

GT  2=THET (2,  1 ,1) - (FLW/SIGMA/CF) **. 25D0 

THET 1=THET (2,1,1) +ERR2 

ERR  2=GT2* (THET (2,1 ,1) -THET1) / (GT2-GT) 

THET  (2,1,1) =TH ET  (2,1,1) -ERR2 

GT=GT2 

IF (DABS (ERR2) . LT.DER2) GOTO  160 
GOTO  100 
160  CONTINUE 
RETURN 

170  CONTINUE 
T2=TO+T 

WRITE  (5,190)  FLW ,T2 ,1 
STOP 
C 

180  FORMAT (3H0**,/,3H  **,115,'  #  ITR  XCD  : 

8  3H  **,F15.3,»  TIME* , /,3  H  **,115 

8  3H  **,F15.6,  '  ERR2* , /, 3H  **) 

190  FORMAT (3H0**,/,3H  **,G15-2,'  FLW  <  0  : 

8  3H  **,F15.3, *  TIME* , /,3H  **,115, 

8  3H  **) 


TSFC  ERR  '  ,/, 

,'  STN* , /, 

TSFC  ERR  * ,/, 

1  STN  * ,  /, 


END 


. 
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SUBROUTINE  PROF 

IMPLICIT  REAL*8 (A-H,0-Z) , INTEGER (I-N) 

DOUBLE  PRECISION  A  ( 50 ) , B  (50 )  , C  (50 ) , D  (50) , F (50)  , G (50)  , 

£  AS  (10)  ,BS  (10)  ,CS  (10)  , 

&  DS  (10)  ,FS (10) ,GS(10)  , 

&  KH  (10)  ,KHH  (50,  10) 

COMMON  /A1/I, J 

8  /A2/TO,TF,DT,T,TM,TL 

£  /A5/SST  (10)  ,  SRT  (10)  ,ZET  ( 1  0)  ,  H  AN  ( 1  0)  ,  DEC,  EOT 

&  /A6/NSTNS,NZPTS,NZPTS1 ,NZPTSA , NPGNO, 

&  NPN0T1,NZSPTS,NZSPT1 

&  /A8/Z  (50)  ,ZS  (10) 

£  /A9/THET (2,50,10) ,TS (2,10,10) 

&  /A10/H22, HI 2,HX, HS22, HS 1 2 , HSX 

&  /A  1  1/LAM, LAMS, SIGMA, CF ,DER2,NITB, NIT2 

&  /A  1 3/KH, KHH 

&  /A14/FLWWV  (10)  ,FSW,FSH,FSHS,FLH,FLW 

&  /A1 5/DZ  A ( 50)  ,DZB  (50) , DZC (50) , C2 (50)  ,C3 (50) 

£  /A 16/DZ AS  (10)  , DZBS (10)  ,DZCS(10)  ,C2S  (10)  ,C3S(10) 

£  /A  17/SO, W,HATO,ALB, DER4 , NITD , NIT4 


100 


T0LD=0.  do 
NIT  4=1 

DO  100  J=  2  ,  NZ PTS  1 

D  ( J)  =C3  (J)  *  (  (THET  (1,  J+1,I)  -THET  (1  ,  J  ,  I)  ) 
£  /DZB  (J)  *KHH(J,I) 

£  -  (THET(1,J,I)  - THE  T ( 1  ,  J-  1,1)  ) 

£  /DZA  (J)  *KHH (J-1 ,1)  ) 

£  *  TH  ET  ( 1  ,  J  ,  I) 

A  ( J)  =KHH  (J,  I)  /DZB(J)  *C2  (J) 

C(J)  =  KHH(J-1,I)/DZA  (J)  *C2  (J) 

B  (J)  =A  (J)  +C  ( J)  +1.  DO 
CONTINUE 


110 


DO  110  J=2 ,NZSPT  1 

D S ( J) =C3S(J) *(  (TS  (1, J+  1,I)-TS (1, J,I)  )  /DZBS (J) 

£  -(TS(1,J,I)-TS(1,J-1,I))/DZAS(J)  ) 

£  +TS ( 1 , J , I ) 

AS  (J)  =  1 .  DO/DZBS  ( J)  *C2  S  (J) 

CS  (J)=1.D0/DZAS(J)  *C2  S  ( J) 

ES  (J) =AS(J)  +CS  (J)  +1.D0 
CONTINUE 


FSW  =  SO*(  DSIN(DEC) *DCOS (ZET(I)  ) 

£  +DCOS  (DEC) *DSIN  (ZET (I) ) * 

£  DCOS  (HAN (I) - W*3600,  DO* (T)  -HATO) 

£  *(1. DO-ALB) 

IF  (  (  (T+TO)  .  LT.  S ST  (I)  )  .AND. 

£  (  (T+  TO)  .  GT.  SRT  (I )  ))  GOTO  120 

IF  (  (  (T  +  TO)  .LT.  (S ST  (I)  +24.  DO)  )  .AND. 

£  (  (T  +  TO)  .GT.  (SRT(I) +24.D0)  )  )  GOTO  120 

FSW=0. DO 


120 


CONTINUE 


(U,r  t*  1  '  *  '  i*  r  •  '  ff 
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CALL  TSFC 

ERR4=THET(2,1,I)  -TOLD 
IF  (  NIT4.LT. 10  ) GOTO  122 
IF  (  NIT4.LT.  50  )  GOTO  125 

THET  (2, 1,1)  =  (  THET  (2, 1,I)+TS (2,1,1)  )/2.D0 

GOTO  125 
122  CONTINUE 

THET  (2,1,1) =THET  (2,1,1)  +  (THET (2,1,1) -TS (2, 1,1) ) *2- DO 
125  CONTINUE 

TS (2, 1,1) =THET (2,1 ,1) 

IF (  DABS (ERR4) -LT. DER4  ) RETURN 
C 

NIT4=NIT4  +  1 

IF  (NII4.LT.  NITD)  GOTO  130 
T  4=TO+T 

WRITE  (5, 18  0)  NIT 4, T4, I , ERR4 
STOP 

130  CONTINUE 
C 

F(NZPTS)  =THET (2, NZPTS,I) 

FS(NZSPTS) =TS  (2,NZSPTS, I) 

G(NZPTS)  =0.  DO 
GS  (NZSPTS)  =0.  DO 
C 

DO  140  J=2 ,NZPTS 1 

JJ  =NZPTS-J  + 1 

E  =  B(JJ)-A(JJ)  *G(  JJ  +  1) 

F(  JJ)  =  (D(JJ)  +  A  (JJ)  *F  (JJ+1)  )  /E 
G  (JJ)  =  C(JJ)  /E 
140  CONTINUE 
C 

DO  150  J=2 ,NZSPT  1 

JJ  =  NZSPTS-J+1 

ES  =  3S (JJ) -AS(JJ)  *GS  (JJ+1) 

FS  ( J  J)  =  (DS  (JJ)  +  AS  (JJ)  *FS(JJ+1)  )  /ES 
GS(JJ)  =  CS(  JJ)  /ES 
150  CONTINUE 
C 

DO  160  J=2,NZPTS1 

THET  (2 ,J, I) =F ( J)  +G  ( J) *T HET  (2, J-1,I) 

160  CONTINUE 
C 

DO  170  J=2 , NZSPT  1 

TS (2, J,I) =FS(J) +GS  (J) *TS  (2,J- 1 ,1) 

170  CONTINUE 
C 

TOL  D= THET (2, 1,1) 

GOTO  120 
C 

180  FORMAT  (3H0**,/,3H  **,115,'  #  ITR  XCD  :  PROF  ERR  ',/, 

&  3H  **,F15.3, '  TIME* ,/,3H  **,115,'  S TN ' , / , 

S  3H  **, F 1 5 . 6 , '  ERR4 ' , /, 3H  **) 


END 
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SUBROUTINE  FVEL 

IMPLICIT  REAL* 8 (A-H,0-Z) , INTEGER (I-N) 

LOGICAL* 1  F  ( 1 ) / 1  *  '/ 

DOUBLE  PRECISION  KN  (1 0)  , KNN (50 , 1 0) , KMOL 
COMMON  /A3/ALAT, DUM1 

B  /A6/NSTNS,NZPTS,NZPTS1 ,IDUM3, NPGNO, IDUM4 (3) 

&  /A8/Z  (50)  ,DUM2  (10) 

&  /A  1 2/KNr  K  NN 

&  /A  1  8/UST AR  ( 1  0)  ,VKC,KMOL 

&  /A19/DUM3 (3)  ,UB 

C 

c 

READ  (3,F)  (USTAR  (I)  ,1  =  1, NSTNS) ,UB 
C 

c 

DO  120  1=1 r NSTNS 

ZT3  =  - 455D 0* USTAR  (I) /2- D0/7. 2  92D-05/DSIN ( AL AT) 

C 

KN  (I) = USTAR  (I) * VKC/2.  DO* Z (NPGNO) 

B  *(  DEXP  (-4.D0*Z (NPGNO)  /ZT3) 

5  +1.D  0/(  1 . D0+ 1 6 .  DO* 

6  (Z (NPGNO) /ZT3) **1. 6D0  )  ) 

&  +  KMCL 


C 


& 

& 

& 

& 

1 1  0 


DO  110  J= 1 ,NZPTS 1 
ZK= (Z (J)  +Z (J+1) )  /2-  DO 
KNN (J, I) =USTAR (I) * VKC/2.  D0*ZK 
*(  DEXP(-4.D0*ZK/ZT3) 

*  1 . DO/  (  1.D0+16. DO* 
(ZK/ZT3) **1.6D0  )  ) 

+  K  MOL 

CONTINUE 


C 

120  CONTINUE 


RETURN 

END 


' 
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SUBROUTINE  DFVTY 

IMPLICIT  REAL*8  (  A-H,0-Z)  ,  INTEGER  (I-N) 

LO  GICA  L* 1  F (1) /'  *•/ 

DOUBLE  PRECISION  KN  ( 1 0)  ,  KH  ( 1  0)  ,  KMOL 
COMMON  /A1/I,J 
8  /A  2/DU  Ml ( 6) 

&  /A8/Z  (50)  ,DUM2  (1  0) 

8  /A9/TH ET(2,50,10)  ,DUM3  (2,10,10) 

&  /A12/KN,  DUM4  (5  0, 10) 

&  /A13/KH,DUM5  (50,  10) 

&  /A18/USTAR  (10)  ,VKC,KMOL 

8  /A19/RI,PHI,BRI,UB 

8  /A2 1 /DUM6 , TSLP , J  SLP ,IDUM2 

C 

TH  ETA V=  (  THET  (1, JSLP- 1,1) +THET ( 1 , JSLP+1 , I)  )/2-D0 
DTH  =THET  (1 , JSLP+1,1) -THET ( 1 , JSLP- 1 , I) 

DZ  =Z  (JSLP+1)  -Z  (JSLP-1) 

BRI  =9. 806/THETA V*DTH/DZ/UB/UB*Z (JSLP)  *Z (JSLP) 
J=2 
C 

CALL  MOFN 
CALL  SLOPE 
C 

KH  (I) =  KN (I) /PHI 
IF  (KH  (I)  .  GT.KMOL) GOTO  100 
KH  (I)  =KMOL 
100  CONTINUE 
C 

RETURN 

END 


SUBROUTINE  SLOPE 

IMPLICIT  REAL* 8 (A-H,0-Z)  , INTEGER (I-N ) 
LOGICAL* 1  F  ( 1 ) / ' *  */ 

COMMON  /A1/I,J 

8  /A2/TO,  DUM4  (2)  ,T,DUN15  (2) 

&  /A  1  8/USTA  R  (  1  0)  ,DUM(2) 

&  /A 19/RI,PHI,BRI,UB 

&  /A2  0/K , ID  UM 

&  /A  21  /DUM2,PHILIM,TSLP ,JSLP,IDUM2 

C 

IF ( T+TO. LT. TSLP) RETURN 
IF  (I.NE.2)  RETURN 
IF  ( J.  GT- JSLP)  RETURN 
READ  (1  ,F)  U 

CF=U*  .4D0/DLOG  (  (.  25D0  +  .  8D  0) /.25D0  ) 

PHI=PHI*UST  AR (I) /CF 
IF (PHI. LT.  PHILIM) RETURN 
PHI=PHILIM 
C 

RETURN 

END 
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SUBROUTINE  DFVTY2 

IMPLICIT  REAL*8  (A-H,0-Z)  , INTEGER  (I-N) 

LO  GIC A  L* 1  F (1) /• *»/ 

DOUBLE  PRECISION  KN  ( 1  0)  ,  KNN  ( 5  0,  1  0)  , 

&  KH(1  0)  ,KHH  (50, 10)  ,KMOL 

COMMON  /A1/I,J 

&  /A2/TO,TF,DT,T,TM,TL 

&  /A6/IDUM2  (2)  ,NZPTS1,IDUM1  (5) 

&  /A8/Z  (50)  ,ZS  (10) 

&  /A9/THET (2,50,10) ,15(2,10,10) 

&  /A  1 2/KN, KNN 

&  /A1 3/KH, KHH 

&  /A  1 8/USTAR (10)  ,VKC,KMOL 

S  /A19/RI,PHI,BRI, UB 


FLAG=0- DO 

KHH  (1 rI)=KNN  (1,1) /PHI 
IF(KHH(1,I)  ,GT,  KMOL)  GOTO  100 
KHH (1,1) =K  MOL 
100  CONTINUE 


DO  130  J=2 , N ZPTS  1 
C 

IF  (  J-LT.10  .OR.  I.EQ.1  )  GOTO  110 
I  F  (  FLAG.  EQ.  1.  DO  )  GOTO  110 
FLAG=1 . DO 
CALL  MCFN 
CALL  SLOPE 
110  CONTINUE 
C 

KHH  (J  ,1)  =KNN  (J,  I)  /PHI 
IF(KHH  (J,I)  -GT.  KMOL)  GOTO  120 
KHH  (J,  I)  =KMOL 
120  CONTINUE 
130  CONTINUE 
C 

RETURN 

END 


' 
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SUBROUTINE  MOFN 

IMPLICIT  REAL*8 (A-H,0-Z) ,  INTEGER (I-N) 

LOGICA L*  1  F  (1)  /'  *'/ 

COMMON  /A  1/1 ,  J 

&  /A2/TO#  DUM4  (2)  ,T#DUM5(2) 

6  /A8/Z  (5  0)  ,  DUM3  (10) 

&  /A19/RI,PHI,BRIr UB 

&  /A 2 1/PHILIM, TS  LP, JSLP, IDUM2 

C 
C 

PHI=. 74D0 

I F  (  (I.  GE.  2,  AND.  J.  LE.  JS  IP)  .AND,  TO+T.GT.TSLP  )  RETURN 
IF  (BRI.EQ. 0. DO)  RETURN 
IF  (BRI. GT. 0. DO) GOTO  150 
C 
C 

N=0 

GZL=BRI 
110  CONTINUE 

X  =  (1 .DO- 15. D0*GZL) **. 25 

PSI=  2 . DO*DLOG (  ( 1 . DO  +  X) /2. DO  ) 

5  +  DLOG  (  (1 .  D0  +  X*X)  /2  .  DO  ) 

6  -2 . DO*  DATAN  (X)  +3.  1 41 592654/2 . DO 
RI  =(  DLOG  (Z  (JSLP) /.0  IDO)  -PSI  )  *X 

R I  =  RI*RI*B  RI 

FNZL=8. 21 4D0*GZL*GZL*GZL-, 5476D0*GZL*GZL-9 .D0*RI*RI*GZL 
&  +R I*RI 

DRI  =(  PSI  — X* (  4. D0/(1.D0+X)  ♦ 1 . DO/ ( 1 . D0+X*X)  )  ) 

&  *15.  DO/4.  D0/(X**3.  DO) 

FNZL  E=-9. D0*RI*DRI*2. DO  +2.D0*RI*DEI 

FNZLP=FNZLP+24. 6  42D0*GZL*GZL- 1. 0952D0*GZL-9.D0*RI*RI 

IF  (FNZLP. EQ.O. DO) GOTO  140 

ERR  =FNZL/FNZLP 

ERRM=ERR/GZL 

GZL  =GZL-ERR 

N  =  N  +  1 

IF  (  DABS  (ERRM)  •  LT.  .  0  5D0)  GOTO  120 
IF(N.GT. 100)  GOTO  130 
GOTO  110 
120  CONTINUE 

PHI  =  .74D0/{  DSQRT  (1.  DO- 9. DO*  GZL)  ) 

RETURN 

C 

130  WRITE  (5,210)  N 
STOP 

140  CONTINUE 

WRITE  (5#220) 

STOP 

C 

C 

150  CONTINUE 

IF  (BRI.  LT.  .  21D0)  GOTO  160 
BRI  =.  2 IDO 
160  CONTINUE 


■ 


, 

n  o 
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A  =22. 09D0*B RI -4.  7D0 

B  =  9.4  DO* DLO G ( 2  (JSLP ) /. 0 1  DO) *BRI - »  74  DO 

C  =BRI*DLOG  (Z  (JSLP)  /.  01  DO)  *  DLOG  (Z  (JSLP)  /.01D0) 

D  =  B*B-  4. DO* AC 

IF (D.LT.O. DO) GOTO  140 

GZL  =  {-B  -DSQRT(D)  )/(2.D0*A) 

PHI  =, 74D0+4. 7D0*GZL 
IF(PHI.LT.PHILIM) RETURN 
PHI  =PHILIM 
RETURN 


210  FORMAT (3H0**,/,3H  **,115,'  #  ITR  XCD  :  MOFN  ERR 

&  3H  **) 

220  FORMAT  (3H0**,/r30H  **  DIVISION  BY  0  :  MOFN  ERR  ,/, 
6  3H  **) 

C 

END 
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IMPLICIT  RE  AL*8 ( A-H, 0-Z) ,  INTEGER (I-N) 

LOGICAL* 1  F  (1 )  /'  ** / 

DIMENSION  Z  (50)  ,  T(50)  ,DZ  (49)  ,  P  (50)  ,  ZS  (6)  ,  TS  (8) 
C 

READ  (5, F)  PB,TAV 

READ  (5,F)  NIDP  ,  (Z(J)  ,T  (J)  ,J=1 ,NIDP) 

READ  (5  ,  F)  NIDPS,  (ZS  (J)  ,TS  (J)  ,J=1,  NIDPS) 

P (NIDP) =PB 

GR  =9 -8  IDO/287.  04D0 

COR  =T (NIDP ) 

NIDP1  =NIDP-1 
C 

DO  100  J= 1 r  NIDP1 
DZ  (J)  =Z  (J)-Z  (J+1) 

100  CONTINUE 
C 

DO  200  J=1,NIDP1 
JJ=NIDP— J 

P  (JJ)  =P  (JJ+  1)  *DEXP  (-GR*  DZ  (JJ)  /TAV) 

200  CONTINUE 
C 

DO  300  J= 1 ,NIDP 
T (J)  =T (J)  * ( (PB/P  (J) )**.  28  6D0) 

300  CONTINUE 
C 

COR=T (NIDP) -COR 
DO  350  J  =  1 f  NIDPS 
TS  (J) =TS (J)  +COR 
350  CONTINUE 
C 

WRITE  (6,4  00)  NIDP,  (Z  (J)  ,T  (J)  ,  J=1,NIDP) 

WRITE  (6,400)  NIDPS,  (ZS(J)  ,TS (J)  ,J=1,NIDPS) 

400  FORMAT  (1 H  , 13 ,  (5  (2F9, 3) ) ) 

C 

STOP 

END 


' 


